Skip to content

Commit

Permalink
Updated srf2stoch to handle mrf files
Browse files Browse the repository at this point in the history
  • Loading branch information
fabiolsilva committed Dec 4, 2024
1 parent 20cd0b3 commit 95c0984
Show file tree
Hide file tree
Showing 3 changed files with 133 additions and 43 deletions.
92 changes: 69 additions & 23 deletions bbp/src/gp/StandRupFormat/srf2stoch.c
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,11 @@ char string[1024];
struct standrupformat srf;
int inbin = 0;

struct srf_apointvalues *apval_ptr;
double dmom, dmu, darea, dslip;
float dd;
int mrf_flag = 0;

float target_dx = -1.0;
float target_dy = -1.0;

Expand Down Expand Up @@ -58,6 +63,10 @@ if(avgstk > -1.0e+14)

read_srf(&srf,infile,inbin);

mrf_flag = 0;
if(strcmp(srf.src_format,"MOMENT") == 0)
mrf_flag = 1;

if(strcmp(outfile,"stdout") == 0)
fpw = stdout;
else
Expand All @@ -73,6 +82,7 @@ sp = NULL;
tr = NULL;
ti = NULL;

apval_ptr = srf.srf_apnts.apntvals;
noff = 0;
for(i=0;i<srf.srf_prect.nseg;i++)
{
Expand Down Expand Up @@ -158,40 +168,76 @@ fprintf(stderr,"id= %d is= %d kp= %d noff= %d\n",id,is,kp,noff);
*/

tt = srf.srf_apnts.apntvals[kp].tinit;
s1 = srf.srf_apnts.apntvals[kp].slip1;
s2 = srf.srf_apnts.apntvals[kp].slip2;

rake = srf.srf_apnts.apntvals[kp].rake;
if(mrf_flag == 1)
{
dmu = (apval_ptr[kp].vs)*(apval_ptr[kp].vs)*(apval_ptr[kp].den);
darea = 1.0*((double)(apval_ptr[kp].area));

cosA = cos(rake*RPERD);
sinA = sin(rake*RPERD);
dmom = sqrt(0.5*apval_ptr[kp].mnn*apval_ptr[kp].mnn
+ 0.5*apval_ptr[kp].mee*apval_ptr[kp].mee
+ 0.5*apval_ptr[kp].mdd*apval_ptr[kp].mdd
+ apval_ptr[kp].mne*apval_ptr[kp].mne
+ apval_ptr[kp].mnd*apval_ptr[kp].mnd
+ apval_ptr[kp].med*apval_ptr[kp].med);

xx = -s2*sinA + s1*cosA;
yy = s2*cosA + s1*sinA;
dslip = dmom/(dmu*darea);
ss = dslip;

ss = sqrt(xx*xx + yy*yy);
/* make a guess on rake based on fault dip */

rr = 90;
if(yy < 0.0)
rr = 270;
dd = apval_ptr[kp].dip;
while(dd < 0.0)
dd = dd + 180.0;
while(dd > 90.0)
dd = 180.0 - dd;

if(xx != 0.0)
{
rr = DPERR*atan(yy/xx);
if(xx < 0.0)
rr = rr + 180;
if(dd <= 45.0)
rr = 90.0;
else
rr = 2.0*(dd-45.0) + 90.0;

tl = (srf.srf_apnts.apntvals[kp].dt)*(srf.srf_apnts.apntvals[kp].ntmr);
}
else
{
s1 = srf.srf_apnts.apntvals[kp].slip1;
s2 = srf.srf_apnts.apntvals[kp].slip2;

rake = srf.srf_apnts.apntvals[kp].rake;

cosA = cos(rake*RPERD);
sinA = sin(rake*RPERD);

while(rr < 0.0)
rr = rr + 360.0;
while(rr > 360.0)
rr = rr - 360.0;
xx = -s2*sinA + s1*cosA;
yy = s2*cosA + s1*sinA;

ss = sqrt(xx*xx + yy*yy);

rr = 90;
if(yy < 0.0)
rr = 270;

if(xx != 0.0)
{
rr = DPERR*atan(yy/xx);
if(xx < 0.0)
rr = rr + 180;
}

while(rr < 0.0)
rr = rr + 360.0;
while(rr > 360.0)
rr = rr - 360.0;

tl = (srf.srf_apnts.apntvals[kp].dt)*(srf.srf_apnts.apntvals[kp].nt1);
if(srf.srf_apnts.apntvals[kp].nt2 > srf.srf_apnts.apntvals[kp].nt1)
tl = (srf.srf_apnts.apntvals[kp].dt)*(srf.srf_apnts.apntvals[kp].nt2);
}

tl = (srf.srf_apnts.apntvals[kp].dt)*(srf.srf_apnts.apntvals[kp].nt1);
if(srf.srf_apnts.apntvals[kp].nt2 > srf.srf_apnts.apntvals[kp].nt1)
tl = (srf.srf_apnts.apntvals[kp].dt)*(srf.srf_apnts.apntvals[kp].nt2);
if(tl < 0.0)
tl = 0.0;

for(k=0;k<nydiv;k++)
{
for(j=0;j<nxdiv;j++)
Expand Down
83 changes: 63 additions & 20 deletions bbp/src/gp/StandRupFormat/srf2stoch_sub.c
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,11 @@ int nx, ny, nxdiv, nydiv, nxsum, nysum;
float shypo, dhypo, elon, elat, fac;
char string[1024];

struct srf_apointvalues *apval_ptr;
double dmom, dmu, darea, dslip;
float dd;
int mrf_flag = 0;

float target_dx = -1.0;
float target_dy = -1.0;

Expand Down Expand Up @@ -49,6 +54,9 @@ if(avgstk > -1.0e+14)
avgstk = avgstk + 360.0;
}

mrf_flag = 0;
if(strcmp(srf->src_format,"MOMENT") == 0)
mrf_flag = 1;

slip = NULL;
trise = NULL;
Expand Down Expand Up @@ -150,6 +158,7 @@ for(i=0;i<srf->srf_prect.nseg;i++)
tr = (float *) check_realloc (tr,nx*ny*sizeof(float));
ti = (float *) check_realloc (ti,nx*ny*sizeof(float));

apval_ptr = srf->srf_apnts.apntvals;
savg = 0.0;
ravg = 0.0;
for(id=0;id<ndip;id++)
Expand All @@ -163,34 +172,68 @@ fprintf(stderr,"id= %d is= %d kp= %d noff= %d\n",id,is,kp,noff);
*/

tt = srf->srf_apnts.apntvals[kp].tinit;
s1 = srf->srf_apnts.apntvals[kp].slip1;
s2 = srf->srf_apnts.apntvals[kp].slip2;

rake = srf->srf_apnts.apntvals[kp].rake;
if(mrf_flag == 1)
{
dmu = (apval_ptr[kp].vs)*(apval_ptr[kp].vs)*(apval_ptr[kp].den);
darea = 1.0*((double)(apval_ptr[kp].area));

dmom = sqrt(0.5*apval_ptr[kp].mnn*apval_ptr[kp].mnn
+ 0.5*apval_ptr[kp].mee*apval_ptr[kp].mee
+ 0.5*apval_ptr[kp].mdd*apval_ptr[kp].mdd
+ apval_ptr[kp].mne*apval_ptr[kp].mne
+ apval_ptr[kp].mnd*apval_ptr[kp].mnd
+ apval_ptr[kp].med*apval_ptr[kp].med);

cosA = cos(rake*RPERD);
sinA = sin(rake*RPERD);
dslip = dmom/(dmu*darea);
ss = dslip;

xx = -s2*sinA + s1*cosA;
yy = s2*cosA + s1*sinA;
/* make a guess on rake based on fault dip */

ss = sqrt(xx*xx + yy*yy);
dd = apval_ptr[kp].dip;
while(dd < 0.0)
dd = dd + 180.0;
while(dd > 90.0)
dd = 180.0 - dd;

rr = 90;
if(yy < 0.0)
rr = 270;
if(dd <= 45.0)
rr = 90.0;
else
rr = 2.0*(dd-45.0) + 90.0;

if(xx != 0.0)
{
rr = DPERR*atan(yy/xx);
if(xx < 0.0)
rr = rr + 180;
tl = (srf->srf_apnts.apntvals[kp].dt)*(srf->srf_apnts.apntvals[kp].ntmr);
}
else
{
s1 = srf->srf_apnts.apntvals[kp].slip1;
s2 = srf->srf_apnts.apntvals[kp].slip2;

rake = srf->srf_apnts.apntvals[kp].rake;

while(rr < 0.0)
rr = rr + 360.0;
while(rr > 360.0)
rr = rr - 360.0;
cosA = cos(rake*RPERD);
sinA = sin(rake*RPERD);

xx = -s2*sinA + s1*cosA;
yy = s2*cosA + s1*sinA;

ss = sqrt(xx*xx + yy*yy);

rr = 90;
if(yy < 0.0)
rr = 270;

if(xx != 0.0)
{
rr = DPERR*atan(yy/xx);
if(xx < 0.0)
rr = rr + 180;
}

while(rr < 0.0)
rr = rr + 360.0;
while(rr > 360.0)
rr = rr - 360.0;
}

tl = (srf->srf_apnts.apntvals[kp].dt)*(srf->srf_apnts.apntvals[kp].nt1);
if(srf->srf_apnts.apntvals[kp].nt2 > srf->srf_apnts.apntvals[kp].nt1)
Expand Down
1 change: 1 addition & 0 deletions bbp/src/gp/StandRupFormat/srf_subs.c
Original file line number Diff line number Diff line change
Expand Up @@ -2039,6 +2039,7 @@ for(i=0;i<mrf[0].srf_apnts.np;i++)
apval_out[i].dt = apval_in[i].dt;
apval_out[i].vp = apval_in[i].vp;
apval_out[i].vs = apval_in[i].vs;
apval_out[i].den = apval_in[i].den;

apval_out[i].stk = apval_in[i].stk;
apval_out[i].dip = apval_in[i].dip;
Expand Down

0 comments on commit 95c0984

Please sign in to comment.