/*-------------------------------------------------------------------- * $Id: psmeca.c,v 1.36 2007/03/24 01:42:07 pwessel Exp $ * * Copyright (c) 1996-2007 by G. Patau * Distributed under the GNU Public Licence * See README file for copying and redistribution conditions. *--------------------------------------------------------------------*/ /* psmeca will read pairs (or ) from inputfile and plot symbols on a map. Focal mechanisms are specified (double couple or moment tensor). PostScript code is written to stdout. Author: Genevieve Patau Date: 7 July 1998 Version: 4 Roots: based on psxy.c Last change 18 February 2000 */ typedef int BOOLEAN; /* BOOLEAN used for logical variables */ typedef int (*PFI) (); /* PFI declares a pointer to a function returning an int */ typedef BOOLEAN (*PFB) (); /* PFB declares a pointer to a function returning a BOOLEAN */ #include #include #include "meca.h" #include "utilmeca.h" int main (int argc, char **argv) { int i, n; /* int i, symbol = 0, n, ix = 0, iy = 1, n_files = 0, fno; */ BOOLEAN outline = FALSE; BOOLEAN plot_zerotrace = FALSE; double xy[2]; double plot_x=0.0, plot_y=0.0, scale = 0.0; double t11 = 1.0, t12 = 0.0, t21 = 0.0, t22 = 1.0; double delaz; char string[BUFSIZ]; char col[15][64]; /* FILE *fp = NULL; struct GMT_PEN pen, lpen, cpen, tr_pen, tz_pen; struct GMT_PEN ppen, tpen; */ int fill[3], efill[3]; /* struct GMT_FILL fill, efill, nofill; struct GMT_FILL pfill, tfill, bfill; */ /* struct nodal_plane NP1, NP2; */ st_me meca; struct MOMENT moment; struct M_TENSOR mt; struct AXIS T, N, P; /* int n_plane = 0, n_plane_old = 0; int default_justify = 2; int justify = default_justify, form = 0; */ int new = 1; /* double default_fontsize = 9.0; double fontsize = default_fontsize; double default_offset = gmtdefs.measure_unit == 0 ? 0.1 : 0.04; double default_pointsize = 0.005;*/ /* double pointsize = default_pointsize, offset = default_offset, angle = 0.0; double fault, depth, depmin = 0.0, depmax = 900.0; */ double size; /*, a_size = GMT_d_NaN; double P_x, P_y, T_x, T_y; char P_sym_type = 0, T_sym_type = 0; int P_sym = 0, T_sym = 0; */ /* fscanf (stdin, "%s %s %s %s %s %s %s %s %s %s %s %[^\n]\n", */ sscanf ("1 1 15 0.38 -7.14 6.76 -0.56 0.32 0.68 23 60 15 081298A", "%s %s %s %s %s %s %s %s %s %s %s %[^\n]\n", col[0], col[1], col[2], col[3], col[4], col[5], col[6], col[7], col[8], col[9], col[10], string); printf("%d\n",argc); for (i=0; i= 360.) str -= 360.; else if(str < 0.) str += 360.; return(str); } /**********************************************************************/ double computed_mw(struct MOMENT moment,double ms) /* Compute mw-magnitude from seismic moment or MS magnitude. */ /* Genevieve Patau from Thorne Lay, Terry C. Wallace Modern Global Seismology Academic Press p. 384 */ { double mw; if(moment.exponent == 0) mw = ms; else mw = (log10(moment.mant) + (double)moment.exponent - 16.1) * 2. / 3.; return(mw); } /*********************************************************************/ double datan2(double y,double x) /* compute arctg in degrees, between -180 et +180. */ /* Genevieve Patau */ { double arctg; double rdeg = 180. / 3.14159265; if(fabs(x) < EPSIL) { if(fabs(y) < EPSIL) { /* fprintf(stderr, "undetermined form 0. / 0."); exit(); */ arctg = 0.; } else arctg = y < 0. ? -90. : 90.; } else if(x < 0.) arctg = y < 0. ? atan(y / x) * rdeg - 180. : atan(y / x) * rdeg + 180.; else arctg = atan(y / x) * rdeg; return(arctg); } /*********************************************************************/ double computed_strike1(struct nodal_plane NP1) /* Compute the strike of the decond nodal plane when are given strike, dip and rake for the first nodal plane with AKI & RICHARD's convention. Angles are in degrees. */ /* Genevieve Patau */ { double str2; double cd1 = cosd(NP1.dip); double temp; double cp2, sp2; double am = NP1.rake == 0 ? 1. : NP1.rake /fabs(NP1.rake); double ss, cs, sr, cr; sincosd (NP1.rake, &sr, &cr); sincosd (NP1.str, &ss, &cs); if(cd1 < EPSIL && fabs(cr) < EPSIL) { /* fprintf(stderr, "\nThe second plane is horizontal;"); fprintf(stderr, "\nStrike is undetermined."); fprintf(stderr, "\nstr2 = NP1.str + 180. is taken to define"); fprintf(stderr, "\nrake in the second plane."); */ str2 = NP1.str + 180.; } else { temp = cr * cs; temp += sr * ss * cd1; sp2 = -am * temp; temp = ss * cr; temp -= sr * cs * cd1; cp2 = am * temp; str2 = datan2(sp2, cp2); str2 = zero_360(str2); } return(str2); } /*********************************************************************/ double computed_dip1(struct nodal_plane NP1) /* Compute second nodal plane dip when are given strike, dip and rake for the first nodal plane with AKI & RICHARD's convention. Angles are in degrees. */ /* Genevieve Patau */ { double am = NP1.rake == 0 ? 1. : NP1.rake / fabs(NP1.rake); double dip2; dip2 = acos(am * sind(NP1.rake) * sind(NP1.dip)) / D2R; return(dip2); } /*********************************************************************/ double computed_rake1(struct nodal_plane NP1) /* Compute rake in the second nodal plane when strike ,dip and rake are given for the first nodal plane with AKI & RICHARD's convention. Angles are in degrees. */ /* Genevieve Patau */ { double computed_strike1(); double computed_dip1(); double rake2, sinrake2; double str2 = computed_strike1(NP1); double dip2 = computed_dip1(NP1); double am = NP1.rake == 0. ? 1. : NP1.rake / fabs(NP1.rake); double sd, cd, ss, cs; sincosd (NP1.dip, &sd, &cd); sincosd (NP1.str - str2, &ss, &cs); if(fabs(dip2 - 90.) < EPSIL) sinrake2 = am * cd; else sinrake2 = -am * sd * cs / cd; rake2 = datan2(sinrake2, -am * sd * ss); return(rake2); } /*********************************************************************/ double computed_dip2(double str1,double dip1,double str2) /* Compute second nodal plane dip when are given strike and dip for the first plane and strike for the second plane. Angles are in degrees. Warning : if dip1 == 90 and cos(str1 - str2) == 0 the second plane dip is undetermined and the only first plane will be plotted. */ /* Genevieve Patau */ { double dip2; double cosdp12 = cosd(str1 - str2); if(fabs(dip1 - 90.) < EPSIL && fabs(cosdp12) < EPSIL) { dip2 = 1000.; /* (only first plane will be plotted) */ } else { dip2 = datan2(cosd(dip1), -sind(dip1) * cosdp12); } return(dip2); } /*********************************************************************/ double computed_rake2(double str1,double dip1,double str2,double dip2,double fault) /* Compute rake in the second nodal plane when strike and dip for first and second nodal plane are given with a double characterizing the fault : +1. inverse fault -1. normal fault. Angles are in degrees. */ /* Genevieve Patau */ { double rake2, sinrake2; double sd, cd, ss, cs; sincosd (str1 - str2, &ss, &cs); sd = sind(dip1); cd = cosd(dip2); if(fabs(dip2 - 90.) < EPSIL) sinrake2 = fault * cd; else sinrake2 = -fault * sd * cs / cd; rake2 = datan2(sinrake2, - fault * sd * ss); return(rake2); } /*********************************************************************/ void define_second_plane(struct nodal_plane NP1,struct nodal_plane *NP2) /* Compute strike, dip, slip for the second nodal plane when are given strike, dip and rake for the first one. */ /* Genevieve Patau */ { NP2->str = computed_strike1(NP1); NP2->dip = computed_dip1(NP1); NP2->rake = computed_rake1(NP1); } /*********************************************************************/ double null_axis_dip(double str1,double dip1,double str2,double dip2) /* compute null axis dip when strike and dip are given for each nodal plane. Angles are in degrees. */ /* Genevieve Patau */ { double den; den = asin(sind(dip1) * sind(dip2) * sind(str1 - str2)) / D2R; if (den < 0.) den = -den; return(den); } /*********************************************************************/ double null_axis_strike(double str1,double dip1,double str2,double dip2) /* Compute null axis strike when strike and dip are given for each nodal plane. Angles are in degrees. */ /* Genevieve Patau */ { double phn, cosphn, sinphn; double sd1, cd1, sd2, cd2, ss1, cs1, ss2, cs2; sincosd (dip1, &sd1, &cd1); sincosd (dip2, &sd2, &cd2); sincosd (str1, &ss1, &cs1); sincosd (str2, &ss2, &cs2); cosphn = sd1 * cs1 * cd2 - sd2 * cs2 * cd1; sinphn = sd1 * ss1 * cd2 - sd2 * ss2 * cd1; if(sind(str1 - str2) < 0.) { cosphn = -cosphn; sinphn = -sinphn; } phn = datan2(sinphn, cosphn); if(phn < 0.) phn += 360.; return(phn); } /*********************************************************************/ double proj_radius(double str1,double dip1,double str) /* Compute the vector radius for a given strike, equal area projection, inferior sphere. Strike and dip of the plane are given. */ /* Genevieve Patau */ { double dip, r; if(fabs(dip1 - 90.) < EPSIL) { /* printf("\nVertical plane : strike is constant."); printf("\nFor ps_mechanism r == 1 for str = str1"); printf("\n else r == 0. is used."); */ r = (str == str1 || str == str1 + 180) ? 1. : 0.; } else { dip = atan(tand(dip1) * sind(str - str1)); r = sqrt(2.) * sin(M_PI_4 - dip / 2.); } return(r); } /*********************************************************************/ /* Compute p-T axis strikes and dips from seismic moment tensor components. Input moment tensor components mrr mtt mff mrt mrf mtf. Angles are in degrees. value, azimuth (en degres), plunge (in degrees) for T, N, P axis. Genevieve Patau, 18 mars 1999 */ #define NP 3 #define NMAT 1 #define NRANSI #define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau);\ a[k][l]=h+s*(g-h*tau); void jacobi(float **a, int n, float d[], float **v, int *nrot) { int j,iq,ip,i; float tresh,theta,tau,t,sm,s,h,g,c,*b,*z; b=vector(1,n); z=vector(1,n); for (ip=1;ip<=n;ip++) { for (iq=1;iq<=n;iq++) v[ip][iq]=0.0; v[ip][ip]=1.0; } for (ip=1;ip<=n;ip++) { b[ip]=d[ip]=a[ip][ip]; z[ip]=0.0; } *nrot=0; for (i=1;i<=50;i++) { sm=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) sm += (float)fabs((double)a[ip][iq]); } if (sm == 0.0) { free_vector(z,1,n); free_vector(b,1,n); return; } if (i < 4) tresh=(float)0.2*sm/((float)(n*n)); else tresh=0.0; for (ip=1;ip<=n-1;ip++) { for (iq=ip+1;iq<=n;iq++) { g=(float)(100.0*fabs(a[ip][iq])); if (i > 4 && (float)(fabs(d[ip])+g) == (float)fabs(d[ip]) && (float)(fabs(d[iq])+g) == (float)fabs(d[iq])) a[ip][iq]=0.0; else if (fabs(a[ip][iq]) > tresh) { h=d[iq]-d[ip]; if ((float)(fabs(h)+g) == (float)fabs(h)) t=(a[ip][iq])/h; else { theta=(float)0.5*h/(a[ip][iq]); t=(float)(1.0/(fabs(theta)+sqrt(1.0+theta*theta))); if (theta < 0.0) t = -t; } c=(float)(1.0/sqrt(1+t*t)); s=t*c; tau=s/((float)1.0+c); h=t*a[ip][iq]; z[ip] -= h; z[iq] += h; d[ip] -= h; d[iq] += h; a[ip][iq]=0.0; for (j=1;j<=ip-1;j++) { ROTATE(a,j,ip,j,iq) } for (j=ip+1;j<=iq-1;j++) { ROTATE(a,ip,j,j,iq) } for (j=iq+1;j<=n;j++) { ROTATE(a,ip,j,iq,j) } for (j=1;j<=n;j++) { ROTATE(v,j,ip,j,iq) } ++(*nrot); } } } for (ip=1;ip<=n;ip++) { b[ip] += z[ip]; d[ip]=b[ip]; z[ip]=0.0; } } nrerror("Too many iterations in routine jacobi"); } #undef ROTATE #undef NRANSI /* (C) Copr. 1986-92 Numerical Recipes Software W".. */ void momten2axe(struct M_TENSOR mt,struct AXIS *T,struct AXIS *N,struct AXIS *P) { /* * routine jacobi from Numerical Recipes is used. * W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery * Numerical Recipes in C * Cambridge University press */ int j,kk,nrot; int jj[3]; float a[3][3]; float *d,*r,**v,**e; float val[3], azi[3], plu[3]; static int num=3; float min,max,mid; float az[3], pl[3]; a[0][0]=(float)mt.f[0]; a[0][1]=(float)mt.f[3]; a[0][2]=(float)mt.f[4]; a[1][0]=(float)mt.f[3]; a[1][1]=(float)mt.f[1]; a[1][2]=(float)mt.f[5]; a[2][0]=(float)mt.f[4]; a[2][1]=(float)mt.f[5]; a[2][2]=(float)mt.f[2]; d=vector(1,NP); r=vector(1,NP); v=matrix(1,NP,1,NP); e=convert_matrix(&a[0][0],1,num,1,num); jacobi(e,num,d,v,&nrot); /* sort eigenvalues */ max = -10000.; for(j=1;j<=num;j++) if(max<=d[j]) {max=d[j];jj[0]=j;} min = 10000.; for(j=1;j<=num;j++) if(min>=d[j]) {min=d[j];jj[2]=j;} mid = 0.; for(j=1;j<=num;j++) if(j!=jj[0]&&j!=jj[2]) jj[1]=j; for (j=1;j<=num;j++) { kk=jj[j-1]; pl[kk]=(float)(asin(- v[1][kk])); az[kk]=(float)(atan2(v[3][kk],- v[2][kk])); if(pl[kk]<=0.) {pl[kk]=-pl[kk]; az[kk]+=(float)(M_PI);} if(az[kk]<0.) az[kk]+=(float)(2.*M_PI); else if(az[kk]>(float)(2.*M_PI)) az[kk]-=(float)(2.*M_PI); pl[kk]*=(float)(180./M_PI); az[kk]*=(float)(180./M_PI); val[j-1] = d[kk]; azi[j-1] = az[kk]; plu[j-1] = pl[kk]; } T->val = (double)val[0]; T->e = mt.expo; T->str = (double)azi[0]; T->dip = (double)plu[0]; N->val = (double)val[1]; N->e = mt.expo; N->str = (double)azi[1]; N->dip = (double)plu[1]; P->val = (double)val[2]; P->e = mt.expo; P->str = (double)azi[2]; P->dip = (double)plu[2]; } /***************************************************************************************/ double ps_tensor(double x0,double y0,double size,struct AXIS T,struct AXIS N,struct AXIS P,int c_rgb[3],int e_rgb[3], int outline, int plot_zerotrace) { int d, b = 1, m; int i, ii, n = 0, j = 1, j2 = 0, j3 = 0; int npoints; /* int lineout = 1; */ int rgb1[3], rgb2[3]; int big_iso = 0; int pvdebug = 0; double a[3], p[3], v[3]; double vi, iso, f; double fir, s2alphan, alphan; double cfi, sfi, can, san; double cpd, spd, cpb, spb, cpm, spm; double cad, sad, cab, sab, cam, sam; double xz, xn, xe; double az = 0., azp = 0., takeoff, r; double azi[3][2]; double x[400], y[400], x2[400], y2[400], x3[400], y3[400]; double xp1[800], yp1[800], xp2[400], yp2[400]; double radius_size; double si, co; FILE *mtout; mtout=fopen("beachball","w"); c_rgb[0]=255; c_rgb[1]=0; c_rgb[2]=0; e_rgb[0]=255; e_rgb[1]=255; e_rgb[2]=255; a[0] = T.str; a[1] = N.str; a[2] = P.str; p[0] = T.dip; p[1] = N.dip; p[2] = P.dip; v[0] = T.val; v[1] = N.val; v[2] = P.val; if(pvdebug){ fprintf(stderr,"\nPV PS_TENSOR A "); fprintf(stderr,"pv size=%f\n",size); fprintf(stderr,"pv c_rgb=%d,%d,%d ",c_rgb[0],c_rgb[1],c_rgb[2]); fprintf(stderr,"pv e_rgb=%d,%d,%d\n",e_rgb[0],e_rgb[1],e_rgb[2]); fprintf(stderr,"pv outline=%d ",outline); fprintf(stderr,"pv plot_zerotrace=%d\n",plot_zerotrace); fprintf(stderr,"pv plot_x/x0=%f ",x0); fprintf(stderr,"pv plot_y/y0=%f\n",y0); fprintf(stderr,"pv a=%f,%f,%f\n",a[0],a[1],a[2]); fprintf(stderr,"pv p=%f,%f,%f\n",p[0],p[1],p[2]); fprintf(stderr,"pv v=%f,%f,%f\n",v[0],v[1],v[2]); } vi = (v[0] + v[1] + v[2]) / 3.; for(i=0; i<=2; i++) v[i] = v[i] - vi; radius_size = size * 0.5; radius_size = 1.0; if(fabs(squared(v[0]) + squared(v[1]) + squared(v[2])) < EPSIL) { if(pvdebug) fprintf(stderr,"PV PS_TENSOR A\n"); /* pure implosion-explosion */ if (vi > 0.) { /* ps_circle(x0, y0, radius_size*2., c_rgb, lineout); */ fprintf(mtout,"0\n"); fprintf(mtout,"0\n"); return(radius_size*2.); } if (vi < 0.) { /* ps_circle(x0, y0, radius_size*2., e_rgb, lineout); */ return(radius_size*2.); fprintf(mtout,"0\n"); fprintf(mtout,"0\n"); } } if(fabs(v[0]) >= fabs(v[2])) { d = 0; m = 2; } else { d = 2; m = 0; } if(plot_zerotrace) vi = 0.; f = - v[1] / v[d]; iso = vi / v[d]; /* Cliff Frohlich, Seismological Research letters, * Vol 7, Number 1, January-February, 1996 * Unless the isotropic parameter lies in the range * between -1 and 1 - f there will be no nodes whatsoever */ if(pvdebug) fprintf(stderr,"PV PS_TENSOR A\n"); if(iso < -1) { if(pvdebug) fprintf(stderr,"pv iso1=%f\n",iso); /* ps_circle(x0, y0, radius_size*2., e_rgb, lineout); */ return(radius_size*2.); } else if(iso > 1-f) { if(pvdebug) fprintf(stderr,"pv iso2=%f\n",iso); /* ps_circle(x0, y0, radius_size*2., c_rgb, lineout); */ return(radius_size*2.); } if(pvdebug) fprintf(stderr,"PV PS_TENSOR B\n"); sincosd (p[d], &spd, &cpd); sincosd (p[b], &spb, &cpb); sincosd (p[m], &spm, &cpm); sincosd (a[d], &sad, &cad); sincosd (a[b], &sab, &cab); sincosd (a[m], &sam, &cam); if(pvdebug) fprintf(stderr,"PV PS_TENSOR C\n"); for(i=0; i<360; i++) { fir = (double) i * D2R; s2alphan = (2. + 2. * iso) / (3. + (1. - 2. * f) * cos(2. * fir)); if(s2alphan > 1.) big_iso++; else { alphan = asin(sqrt(s2alphan)); sincos (fir, &sfi, &cfi); sincos (alphan, &san, &can); xz = can * spd + san * sfi * spb + san * cfi * spm; xn = can * cpd * cad + san * sfi * cpb * cab + san * cfi * cpm * cam; xe = can * cpd * sad + san * sfi * cpb * sab + san * cfi * cpm * sam; if(fabs(xn) < EPSIL && fabs(xe) < EPSIL) { takeoff = 0.; az = 0.; } else { az = atan2(xe, xn); if(az < 0.) az += M_PI * 2.; takeoff = acos(xz / sqrt(xz * xz + xn * xn + xe * xe)); } if(takeoff > M_PI_2) { takeoff = M_PI - takeoff; az += M_PI; if(az > M_PI * 2.) az -= M_PI * 2.; } r = M_SQRT2 * sin(takeoff / 2.); sincos (az, &si, &co); if(i == 0) { azi[i][0] = az; x[i] = x0 + radius_size * r * si; y[i] = y0 + radius_size * r * co; azp = az; } else { if(fabs(fabs(az - azp) - M_PI) < D2R * 10.) { azi[n][1] = azp; azi[++n][0] = az; } if(fabs(fabs(az -azp) - M_PI * 2.) < D2R * 2.) { if(azp < az) azi[n][0] += M_PI * 2.; else azi[n][0] -= M_PI * 2.; } switch (n) { case 0 : x[j] = x0 + radius_size * r * si; y[j] = y0 + radius_size * r * co; j++; break; case 1 : x2[j2] = x0 + radius_size * r * si; y2[j2] = y0 + radius_size * r * co; j2++; break; case 2 : x3[j3] = x0 + radius_size * r * si; y3[j3] = y0 + radius_size * r * co; j3++; break; } azp = az; } } } azi[n][1] = az; if(v[1] < 0.) for(i=0;i<=2;i++) {rgb1[i] = c_rgb[i]; rgb2[i] = e_rgb[i];} else for(i=0;i<=2;i++) {rgb1[i] = e_rgb[i]; rgb2[i] = c_rgb[i];} if(pvdebug) fprintf(stderr,"PV PS_TENSOR, n=%d\n",n); /* pv makes circle needed by ps: ps_circle(x0, y0, radius_size*2., rgb2, lineout); */ switch(n) { case 0 : if(pvdebug) fprintf(stderr,"PV PS_TENSOR case 0\n"); for(i=0; i<360; i++) { xp1[i] = x[i]; yp1[i] = y[i]; } npoints = i; /* ps_polygon(xp1, yp1, npoints, rgb1, outline); */ fprintf(mtout,"%d\n",npoints); for(i=0; i M_PI) azi[0][0] -= M_PI * 2.; else if(azi[0][1] - azi[0][0] > M_PI) azi[0][0] += M_PI * 2.; if(azi[0][0] < azi[0][1]) for(az = azi[0][1] - D2R; az > azi[0][0]; az -= D2R) { sincos (az, &si, &co); xp1[i] = x0 + radius_size * si; yp1[i++] = y0 + radius_size * co; } else for(az = azi[0][1] + D2R; az < azi[0][0]; az += D2R) { sincos (az, &si, &co); xp1[i] = x0 + radius_size * si; yp1[i++] = y0 + radius_size * co; } npoints = i; /* ps_polygon(xp1, yp1, npoints, rgb1, outline); */ fprintf(mtout,"%d\n",npoints); for(i=0; i M_PI) azi[1][0] -= M_PI * 2.; else if(azi[1][1] - azi[1][0] > M_PI) azi[1][0] += M_PI * 2.; if(azi[1][0] < azi[1][1]) for(az = azi[1][1] - D2R; az > azi[1][0]; az -= D2R) { sincos (az, &si, &co); xp2[i] = x0 + radius_size * si; yp2[i++] = y0 + radius_size * co; } else for(az = azi[1][1] + D2R; az < azi[1][0]; az += D2R) { sincos (az, &si, &co); xp2[i] = x0 + radius_size * si; yp2[i++] = y0 + radius_size * co; } npoints = i; /* ps_polygon(xp2, yp2, npoints, rgb1, outline); */ fprintf(mtout,"%d\n",npoints); for(i=0; i=0; ii--) { xp1[i] = x2[ii]; yp1[i++] = y2[ii]; } npoints = i; /* ps_polygon(xp1, yp1, npoints, rgb1, outline); */ break; } if(azi[2][0] - azi[0][1] > M_PI) azi[2][0] -= M_PI * 2.; else if(azi[0][1] - azi[2][0] > M_PI) azi[2][0] += M_PI * 2.; if(azi[2][0] < azi[0][1]) for(az = azi[0][1] - D2R; az > azi[2][0]; az -= D2R) { sincos (az, &si, &co); xp1[i] = x0+ radius_size * si; yp1[i++] = y0+ radius_size * co; } else for(az = azi[0][1] + D2R; az < azi[2][0]; az += D2R) { sincos (az, &si, &co); xp1[i] = x0+ radius_size * si; yp1[i++] = y0+ radius_size * co; } npoints = i; /* fprintf(stderr,"pv npoints=%d\n",i); */ fprintf(mtout,"%d\n",npoints); for(i=0; i M_PI) azi[1][0] -= M_PI * 2.; else if(azi[1][1] - azi[1][0] > M_PI) azi[1][0] += M_PI * 2.; if(azi[1][0] < azi[1][1]) for(az = azi[1][1] - D2R; az > azi[1][0]; az -= D2R) { sincos (az, &si, &co); xp2[i] = x0+ radius_size * si; yp2[i++] = y0+ radius_size * co; } else for(az = azi[1][1] + D2R; az < azi[1][0]; az += D2R) { sincos (az, &si, &co); xp2[i] = x0+ radius_size * si; yp2[i++] = y0+ radius_size * co; } npoints = i; fprintf(mtout,"%d\n",npoints); for(i=0; i M_PI_2) { d1 = M_PI - d1; p1 += M_PI; if(p1 > PII) p1 -= PII; } if(p1 < 0.) p1 += PII; amz = sdt - sdp; amx = spt - spp; amy = cpt - cpp; d2 = atan2(sqrt(amx*amx + amy*amy), amz); p2 = atan2(amy, -amx); if(d2 > M_PI_2) { d2 = M_PI - d2; p2 += M_PI; if(p2 > PII) p2 -= PII; } if(p2 < 0.) p2 += PII; NP1->dip = d1 / D2R; NP1->str = p1 / D2R; NP2->dip = d2 / D2R; NP2->str = p2 / D2R; im = 1; if(dp > dt) im = -1; NP1->rake = computed_rake2(NP2->str,NP2->dip,NP1->str,NP1->dip,im); NP2->rake = computed_rake2(NP1->str,NP1->dip,NP2->str,NP2->dip,im); } /*********************************************************************/ void dc_to_axe(st_me meca,struct AXIS *T,struct AXIS *N,struct AXIS *P) /* From FORTRAN routines of Anne Deschamps : compute azimuth and plungement of P-T axis from nodal plane strikes, dips and rakes. */ { int im = 0; int pure_strike_slip = 0; double cd1, sd1, cd2, sd2; double cp1, sp1, cp2, sp2; double amz, amx, amy, dx, px, dy, py; if(fabs(sind(meca.NP1.rake)) > EPSIL) im = (int) (meca.NP1.rake / fabs(meca.NP1.rake)); else if(fabs(sind(meca.NP2.rake)) > EPSIL) im = (int) (meca.NP2.rake / fabs(meca.NP2.rake)); else pure_strike_slip = 1; if(pure_strike_slip) { if(cosd(meca.NP1.rake) < 0.) { P->str = zero_360(meca.NP1.str + 45.); T->str = zero_360(meca.NP1.str - 45.); } else { P->str = zero_360(meca.NP1.str - 45.); T->str = zero_360(meca.NP1.str + 45.); } P->dip = 0.; T->dip = 0.; } else { cd1 = cosd(meca.NP1.dip) * M_SQRT2; sd1 = sind(meca.NP1.dip) * M_SQRT2; cd2 = cosd(meca.NP2.dip) * M_SQRT2; sd2 = sind(meca.NP2.dip) * M_SQRT2; cp1 = - cosd(meca.NP1.str) * sd1; sp1 = sind(meca.NP1.str) * sd1; cp2 = - cosd(meca.NP2.str) * sd2; sp2 = sind(meca.NP2.str) * sd2; amz = - (cd1 + cd2); amx = - (sp1 + sp2); amy = cp1 + cp2; dx = atan2(sqrt(amx * amx + amy * amy), amz) - M_PI_2; px = atan2(amy, - amx); if(px < 0.) px += TWO_PI; amz = cd1 - cd2; amx = sp1 - sp2; amy = - cp1 + cp2; dy = atan2(sqrt(amx * amx + amy * amy), - fabs(amz)) - M_PI_2; py = atan2(amy, - amx); if(amz > 0.) py -= M_PI; if(py < 0.) py += TWO_PI; if(im == 1) { P->dip = dy; P->str = py; T->dip = dx; T->str = px; } else { P->dip = dx; P->str = px; T->dip = dy; T->str = py; } } T->str /= D2R; T->dip /= D2R; P->str /= D2R; P->dip /= D2R; N->str = null_axis_strike(T->str, T->dip, P->str, P->dip); N->dip = null_axis_dip(T->str, T->dip, P->str, P->dip); } /*********************************************************************/ void axis2xy(double x0,double y0,double size,double pp,double dp,double pt,double dt,double *xp,double *yp,double *xt,double *yt) /* angles are in degrees */ { double radius; double spp, cpp, spt, cpt; sincosd (pp, &spp, &cpp); sincosd (pt, &spt, &cpt); size *= 0.5; radius = sqrt(1. - sind(dp)); if(radius >= 0.97) radius = 0.97; *xp = radius * spp * size + x0; *yp = radius * cpp * size + y0; radius = sqrt(1. - sind(dt)); if(radius >= 0.97) radius = 0.97; *xt = radius * spt * size + x0; *yt = radius * cpt * size + y0; } /*****************************************************************************/ /* $Id: nrutil.c,v 1.3 2007/02/02 16:23:20 pwessel Exp $ * Public Domain NR stuff. */ #define _POSIX_SOURCE 1 #include #include #include #define NR_END 1 #define FREE_ARG char* void nrerror(char error_text[]) /* Numerical Recipes standard error handler */ { fprintf(stderr,"Numerical Recipes run-time error...\n"); fprintf(stderr,"%s\n",error_text); fprintf(stderr,"...now exiting to system...\n"); exit(1); } float *vector(long nl, long nh) /* allocate a float vector with subscript range v[nl..nh] */ { float *v; v=(float *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(float))); if (!v) nrerror("allocation failure in vector()"); return v-nl+NR_END; } int *ivector(long nl, long nh) /* allocate an int vector with subscript range v[nl..nh] */ { int *v; v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int))); if (!v) nrerror("allocation failure in ivector()"); return v-nl+NR_END; } unsigned char *cvector(long nl, long nh) /* allocate an unsigned char vector with subscript range v[nl..nh] */ { unsigned char *v; v=(unsigned char *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(unsigned char))); if (!v) nrerror("allocation failure in cvector()"); return v-nl+NR_END; } unsigned long *lvector(long nl, long nh) /* allocate an unsigned long vector with subscript range v[nl..nh] */ { unsigned long *v; v=(unsigned long *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(long))); if (!v) nrerror("allocation failure in lvector()"); return v-nl+NR_END; } double *dvector(long nl, long nh) /* allocate a double vector with subscript range v[nl..nh] */ { double *v; v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double))); if (!v) nrerror("allocation failure in dvector()"); return v-nl+NR_END; } float **matrix(long nrl, long nrh, long ncl, long nch) /* allocate a float matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; float **m; /* allocate pointers to rows */ m=(float **) malloc((size_t)((nrow+NR_END)*sizeof(float*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(float *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(float))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } double **dmatrix(long nrl, long nrh, long ncl, long nch) /* allocate a double matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; double **m; /* allocate pointers to rows */ m=(double **) malloc((size_t)((nrow+NR_END)*sizeof(double*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(double *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(double))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } int **imatrix(long nrl, long nrh, long ncl, long nch) /* allocate a int matrix with subscript range m[nrl..nrh][ncl..nch] */ { long i, nrow=nrh-nrl+1,ncol=nch-ncl+1; int **m; /* allocate pointers to rows */ m=(int **) malloc((size_t)((nrow+NR_END)*sizeof(int*))); if (!m) nrerror("allocation failure 1 in matrix()"); m += NR_END; m -= nrl; /* allocate rows and set pointers to them */ m[nrl]=(int *) malloc((size_t)((nrow*ncol+NR_END)*sizeof(int))); if (!m[nrl]) nrerror("allocation failure 2 in matrix()"); m[nrl] += NR_END; m[nrl] -= ncl; for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol; /* return pointer to array of pointers to rows */ return m; } float **submatrix(float **a, long oldrl, long oldrh, long oldcl, long oldch, long newrl, long newcl) /* point a submatrix [newrl..][newcl..] to a[oldrl..oldrh][oldcl..oldch] */ { long i,j,nrow=oldrh-oldrl+1,ncol=oldcl-newcl; float **m; /* allocate array of pointers to rows */ m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); if (!m) nrerror("allocation failure in submatrix()"); m += NR_END; m -= newrl; /* set pointers to rows */ for(i=oldrl,j=newrl;i<=oldrh;i++,j++) m[j]=a[i]+ncol; /* return pointer to array of pointers to rows */ return m; } float **convert_matrix(float *a, long nrl, long nrh, long ncl, long nch) /* allocate a float matrix m[nrl..nrh][ncl..nch] that points to the matrix declared in the standard C manner as a[nrow][ncol], where nrow=nrh-nrl+1 and ncol=nch-ncl+1. The routine should be called with the address &a[0][0] as the first argument. */ { long i,j,nrow=nrh-nrl+1,ncol=nch-ncl+1; float **m; /* allocate pointers to rows */ m=(float **) malloc((size_t) ((nrow+NR_END)*sizeof(float*))); if (!m) nrerror("allocation failure in convert_matrix()"); m += NR_END; m -= nrl; /* set pointers to rows */ m[nrl]=a-ncl; for(i=1,j=nrl+1;i