#include #include #define MAXPHASE 720 /* maximum number of velest phases per event */ #define MAXSTAS 2000 /* maximum number of stations in table */ #define MINPHASE 4 /* minimum number of phases to output per event */ /* a program to convert hypoinverse arc files to convert format for 3d inversion */ /* or velest format */ /* major rewrite to build a table of P and S readings for each event before outputing them */ main() { FILE *fpin,*fpout,*fpout3,*fpfinal; char namein[100],nameout[100],nameout2[100],nameout3[100],namefinal[100]; char string[20]; char line[200]; char output[200]; char command[300],command2[300]; char summout[300]; int informat; /* input format */ int outformat; /* output format */ int count; int nphase; /* number of phases in an event */ int minphase; /* minimum number of phases to output */ int numevents; /* count of events */ int rtp,swaves,fours; /* if 1 use rtp and/or s waves and/or four weights */ int phoriz,svert; /* use P-waves from horizontals or s waves from verticals */ int codamag; /* if one prefer coda mags */ int latd,lond; float latm,lonm,depth,mag; float lat,lon; int yre,mthe,daye,doye,hre,mine; /* date and time of epicenter */ float sece; /* seconds of epicenter */ int yrp,mthp,dayp,doyp,hrp,minp; /* date and time of pick */ int century; /* what century is it */ float secp,secs; /* seconds of p and s picks */ float seca; /* adjusted arrival time */ char station[10],pcomm[10],scom[10],source,comp; int pweight,sweight; int inevent; /* 1 while in an event */ int i,j,k; float x,y,z; int doy(); float adjtime(); char network[5]; /* network identifier */ char newstas[100][6]; int nnew; int countall; char namenew[100]; FILE *fpnew; int goodsta; int whichsta; /* a table for each event to store the good picks */ struct pickstable { char station[10]; int pweight; int sweight; float secp; float secs; float secpa; float secsa; char pcomm[10]; char scomm[10]; } picks[MAXSTAS]; int nstas; /* the number of stations for this event */ numevents=0; /* get input file and make sure it is openable */ i=1; while(i){ sprintf(namein," "); getstr("enter input hypoinverse arc file",namein,1); fpin=fopen(namein,"r"); if(fpin==NULL)fprintf(stderr,"unable to open %s\n",namein); else i=0; } /* is input new or old style ncsn phase file */ sprintf(string,"y"); getstr("Input in y2k, new, or old format (y, n or o)",string,0); if(string[0]=='n')informat=1; if(string[0]=='y')informat=2; else informat=0; /* get output file and make sure it is openable */ i=1; while(i){ sprintf(namefinal," "); getstr("enter output convert format file",namefinal,1); fpout=fopen(namefinal,"w"); if(fpout==NULL)fprintf(stderr,"unable to open %s\n",namefinal); else i=0; } /* should output be velest or 3-d */ sprintf(string,"3"); getstr("Output in 3-d or velest format (3 or v)",string,0); if(string[0]=='v')outformat=0; else outformat=1; /* get station output file and make sure it is openable */ i=1; while(i){ sprintf(nameout2," "); getstr("enter station list output file",nameout2,1); fpout3=fopen(nameout2,"w"); if(fpout3==NULL)fprintf(stderr,"unable to open %s\n",nameout2); else i=0; } fclose(fpout3); /* get temporary station output file and make sure it is openable */ i=1; while(i){ sprintf(nameout3,"%s.temp",nameout2); getstr("enter temporary station list output file",nameout3,0); fpout3=fopen(nameout3,"w"); if(fpout3==NULL)fprintf(stderr,"unable to open %s\n",nameout3); else i=0; } sprintf(command,"cat %s >> %s",nameout,namefinal); sprintf(command2,"cat %s",nameout); /* find minimum number of phases */ minphase=MINPHASE; getint("Enter minimum number of phases to use",&minphase,0); /* get new stations read in */ i=1; while(i){ sprintf(namenew,"all"); getstr("enter all to count phases from all stations or the name of a station list file that do count",namenew,1); if(strcmp("all",namenew)==0){ countall=1; i=0; } else{ countall=0; fpnew=fopen(namenew,"r"); if(fpnew==NULL){ fprintf(stderr,"unable to open %s\n",namenew); continue; } else i=0; } if(countall==0){ nnew=0; while(fscanf(fpnew,"%s",&newstas[nnew][0])!=EOF)nnew++; fclose(fpnew); } } /* find out if we want rtp readings */ sprintf(string,"no"); getstr("Use RTP readings? (yes or no)",string,0); if(string[0]=='y')rtp=1; else rtp=0; /* find out if we want S readings */ sprintf(string,"no"); getstr("Use S waves? (yes or no)",string,0); if(string[0]=='y')swaves=1; else swaves=0; /* do we want S waves from verticals */ if(swaves){ sprintf(string,"no"); getstr("Use S waves from verticals? (yes or no)",string,0); if(string[0]=='y')svert=1; else svert=0; } else svert=0; /* find out if we want P readings from horizontals */ sprintf(string,"no"); getstr("Use P waves from horizontals? (yes or no)",string,0); if(string[0]=='y')phoriz=1; else phoriz=0; /* find out if we want to carry 4 and 9 weights */ sprintf(string,"no"); getstr("Carry weights 4 and above? (yes or no)",string,0); if(string[0]=='y')fours=1; else fours=0; /* find out if we want to carry amp or coda mags */ sprintf(string,"coda"); getstr("Prefer coda or amplitude mags? (coda or amp)",string,0); if(string[0]=='a')codamag=0; else codamag=1; /* sample input file old style 8910180028452536 5642121 4294 1483 0 42 77 13 29 9985 125247 3 5537SJB 32 1 55 12440 0275 0 22LOMT22 2 49 BSRM P 4E891018 028 99994799 0 5639ES 4 -88 0 0 0 35110105 0 3615033 0 0 0 2 BSRM P 4N891018 028 99994799 0 5712ES 3 -15 41 0 0 35110105 0 3015031 0 0 20 2 BSRM P 4Z891018 028 99994799 0 5641ES 4 -85 0 0 0 35110105 0 3015031 0 0 0 2 145642 new style 9701022236092639 98123 369 40438 37 55 5 10 4387 68233 2 1834MAA 16 0 18 6889175999 20 20PTA 222D119X 30124970D342999Z374 20 GHL IPD0V97 1 22236 1020 -8124 0 0 0 0 -20 0 4512600 0 81 5333 0 627 0 2D VHZNC GHL YPD5V97 1 22236 1022 -6 0 0 0 0 0 -20 0 4512605 0 136 5338 0 0 0 PD VHZNC GCW YPU8V97 1 22236 1066 -94 0 0 0 0 0 -29 0 128 9805 0 12535538 0 0 0 PD VHZNC GCW IPD0V97 1 22236 1157 -3124 0 0 0 0 -29 0 128 9800 0 7935534 0 214 0 2D VHZNC y2k style 199309030958211436 5774121 3529 251 00 34 33 1 912686 40 51 0 18160SAR 15 0 18 40 38 0 313 0 12LOMT22 D 38X 30027461D159 313 HCR NC VVHZ XPD6199309030958 2163 -15 0 0 0 0 00 0 0 -3 0 715805 0 23147220 00 0 0RD HCR NC VVHZ IPU0199309030958 2179 1138 0 0 0 00 0 0 -3 0 715800 0 21147210 00 821 02D HCB NC VVHZ IPU0199309030958 2290 0138 0 0 0 00 0 0 1 0 73 8700 0 18242170 00 162 02D HCB NC VVHZ XPU5199309030958 2291 1 0 0 0 0 00 0 0 1 0 73 8705 0 19242180 00 0 0RD HOR NC VVHZ IPU0199309030958 2296 -17138 0 0 0 00 0 0 10 0 81 8700 0 20126150 00 155 02D HOR NC VVHZ XPU5199309030958 2298 -15 0 0 0 0 00 0 0 10 0 81 8705 0 21126160 00 0 0RD HGW NC VVHZ IPD0199309030958 2297 0138 0 0 0 00 0 0 -12 0 84 8700 0 12317160 00 128 02D HGW NC VVHZ XPD5199309030958 2299 2 0 0 0 0 00 0 0 -12 0 84 8705 0 13317160 00 0 0RD HAZ NC VVHZ IPD0199309030958 2317 32 0 0 0 0 00 0 0 -29 0 86 9000 0 12183180 00 0 02D */ /* sample output file for 3-D (outformat=1) 84 530 1623 9.05 35 47.20 120 19.63 3.23 1.20 PGHMIPD0 1.50PMRMIPU0 2.12PMCMIPD0 2.17PHAMIPU0 2.21PSRMIPD0 2.26PAGMIPU0 2.34 PPFMIPD0 3.04PTRMEP 2 3.33PWKMIPU0 3.82PHGMIPU0 3.82PCAMIPD0 3.99PSTMEPD2 4.67 PSHMIPD0 4.75PMGMEPU2 7.66PPTMEPU2 9.31 0 alternate data line for velest (outformat=0) 840210 23 6 50.92 37.1653N 121.5500W 5.62 0.00 73 0.0 0.03 1.0 1.0 GHPFP0 1.41GHPFS4 2.24EDPFP0 1.98EDPFS4 3.38JNPFP0 3.00JNPFS4 5.27 */ inevent=0; while(fgets(line,200,fpin)!=NULL){ if(line[0]=='$')continue; /* skip all shadow cards */ if(!inevent){ /* read the summary line into epicenter info */ if(informat<2){ /* non-y2k formats */ string[0]=line[0]; string[1]=line[1]; string[2]='\0'; sscanf(string,"%d",&yre); yre+=1900; string[0]=line[2]; string[1]=line[3]; string[2]='\0'; sscanf(string,"%d",&mthe); string[0]=line[4]; string[1]=line[5]; string[2]='\0'; sscanf(string,"%d",&daye); string[0]=line[6]; string[1]=line[7]; string[2]='\0'; sscanf(string,"%d",&hre); string[0]=line[8]; string[1]=line[9]; string[2]='\0'; sscanf(string,"%d",&mine); string[0]=line[10]; string[1]=line[11]; string[2]=line[12]; string[3]=line[13]; string[4]='\0'; sscanf(string,"%f",&sece); string[0]=line[14]; string[1]=line[15]; string[2]='\0'; sscanf(string,"%d",&latd); string[0]=line[17]; string[1]=line[18]; string[2]=line[19]; string[3]=line[20]; string[4]='\0'; sscanf(string,"%f",&latm); string[0]=line[21]; string[1]=line[22]; string[2]=line[23]; string[3]='\0'; sscanf(string,"%d",&lond); string[0]=line[25]; string[1]=line[26]; string[2]=line[27]; string[3]=line[28]; string[4]='\0'; sscanf(string,"%f",&lonm); string[0]=line[29]; string[1]=line[30]; string[2]=line[31]; string[3]=line[32]; string[4]=line[33]; string[5]='\0'; sscanf(string,"%f",&depth); if(codamag){ /* sscanf(&line[67],"%2f",&mag);*/ string[0]=line[67]; string[1]=line[68]; string[2]='\0'; sscanf(string,"%f",&mag); } if(!codamag || mag==0){ /* sscanf(&line[34],"%2f",&mag);*/ string[0]=line[34]; string[1]=line[35]; string[2]='\0'; sscanf(string,"%f",&mag); } if(mag==0 && !codamag){ /* sscanf(&line[67],"%2f",&mag);*/ string[0]=line[67]; string[1]=line[68]; string[2]='\0'; sscanf(string,"%f",&mag); } sece/=100; latm/=100; lonm/=100; lat=latd + latm/60.; lon=lond + lonm/60.; depth/=100; mag/=10; doye=doy(yre,mthe,daye); } else { /* y2k format */ string[0]=line[0]; string[1]=line[1]; string[2]=line[2]; string[3]=line[3]; string[4]='\0'; sscanf(string,"%d",&yre); string[0]=line[4]; string[1]=line[5]; string[2]='\0'; sscanf(string,"%d",&mthe); string[0]=line[6]; string[1]=line[7]; string[2]='\0'; sscanf(string,"%d",&daye); string[0]=line[8]; string[1]=line[9]; string[2]='\0'; sscanf(string,"%d",&hre); string[0]=line[10]; string[1]=line[11]; string[2]='\0'; sscanf(string,"%d",&mine); string[0]=line[12]; string[1]=line[13]; string[2]=line[14]; string[3]=line[15]; string[4]='\0'; sscanf(string,"%f",&sece); string[0]=line[16]; string[1]=line[17]; string[2]='\0'; sscanf(string,"%d",&latd); string[0]=line[19]; string[1]=line[20]; string[2]=line[21]; string[3]=line[22]; string[4]='\0'; sscanf(string,"%f",&latm); string[0]=line[23]; string[1]=line[24]; string[2]=line[25]; string[3]='\0'; sscanf(string,"%d",&lond); string[0]=line[27]; string[1]=line[28]; string[2]=line[29]; string[3]=line[30]; string[4]='\0'; sscanf(string,"%f",&lonm); string[0]=line[31]; string[1]=line[32]; string[2]=line[33]; string[3]=line[34]; string[4]=line[33]; string[5]='\0'; sscanf(string,"%f",&depth); if(codamag){ string[0]=line[70]; string[1]=line[71]; string[2]=line[72]; string[3]='\0'; sscanf(string,"%f",&mag); } if(!codamag || mag==0){ string[0]=line[36]; string[1]=line[37]; string[2]=line[38]; string[3]='\0'; sscanf(string,"%f",&mag); } if(mag==0 && !codamag){ string[0]=line[70]; string[1]=line[71]; string[2]=line[72]; string[3]='\0'; sscanf(string,"%f",&mag); } sece/=100; latm/=100; lonm/=100; lat=latd + latm/60.; lon=lond + lonm/60.; depth/=100; mag/=100; doye=doy(yre,mthe,daye); } /* get rid of century to print out year */ century=yre/100; century*=100; if(outformat)sprintf(summout,"%2d%2d%2d %2d%2d %5.2f %2d %5.2f %3d %5.2f %6.2f %5.2f\n", yre-century,mthe,daye,hre,mine,sece,latd,latm,lond,lonm,depth,mag); else sprintf(summout,"%2d%2d%2d %2d%2d %5.2f %7.4fN %8.4fW%7.2f%7.2f 0 \n", yre-century,mthe,daye,hre,mine,sece,lat,lon,depth,mag); /* velest format */ inevent=1; count=0; nphase=0; nstas=0; strcpy(&(picks[0].station[0])," "); continue; } if(inevent){ if(line[0]==' ' || line[0]=='\n' || line[0]=='\0' || line[0]==' '){ /* time to output event */ /* first see if there are enough phases to justify outputting event */ nphase=0; for(i=0;iminphase){ /* actually output the event */ numevents++; /* keep count of output events */ fputs(summout,fpout); /* the summary line */ count=0; for(i=0;i=1){ /* convert new station names to a 4 character code */ if(informat==1){ network[0]=line[98]; network[1]=line[99]; network[2]='\0'; } else { network[0]=line[5]; network[1]=line[6]; network[2]='\0'; } if(strcmp(network,"BK")==0)station[3]='B'; else if(strcmp(network,"CI")==0)station[3]='C'; else if(strcmp(network,"TS")==0)station[3]='C'; else if(strcmp(network,"SN")==0)station[3]='D'; else if(strcmp(network,"PX")==0)station[3]='F'; else if(strcmp(network,"NC")==0)station[3]='M'; else if(strcmp(network,"LL")==0)station[3]='L'; else if(strcmp(network,"UW")==0)station[3]='O'; else if(strcmp(network,"NN")==0)station[3]='R'; else if(strcmp(network,"SC")==0)station[3]='S'; else if(strcmp(network,"UU")==0)station[3]='U'; else if(strcmp(network,"WR")==0)station[3]='W'; else if(strcmp(network,"MX")==0)station[3]='X'; else if(strcmp(network,"PG")==0)station[3]='Y'; else if(strcmp(network,"AZ")==0)station[3]=' '; else station[3]=' '; } if(informat<2){ /* old and new phase lines */ pcomm[0]=line[4]; pcomm[1]=line[5]; pcomm[2]=line[6]; pcomm[3]='\0'; string[0]=line[7]; string[1]='\0'; sscanf(string,"%d",&pweight); if(informat==1)comp=line[97]; else comp=line[8]; string[0]=line[9]; string[1]=line[10]; string[2]='\0'; sscanf(string,"%d",&yrp); yrp+=1900; string[0]=line[11]; string[1]=line[12]; string[2]='\0'; sscanf(string,"%d",&mthp); string[0]=line[13]; string[1]=line[14]; string[2]='\0'; sscanf(string,"%d",&dayp); string[0]=line[15]; string[1]=line[16]; string[2]='\0'; sscanf(string,"%d",&hrp); string[0]=line[17]; string[1]=line[18]; string[2]='\0'; sscanf(string,"%d",&minp); for(i=0;i<5;++i)string[i]=line[i+19]; string[5]='\0'; sscanf(string,"%f",&secp); for(i=0;i<5;++i)string[i]=line[i+31]; string[5]='\0'; sscanf(string,"%f",&secs); for(i=0;i<3;++i)scom[i]=line[36+i]; scom[3]='\0'; string[0]=line[39]; string[1]='\0'; sscanf(string,"%d",&sweight); source=line[91]; secp/=100; secs/=100; doyp=doy(yrp,mthp,dayp); } else { /* y2k phase lines */ pcomm[0]=line[13]; pcomm[1]=line[14]; pcomm[2]=line[15]; pcomm[3]='\0'; string[0]=line[16]; string[1]='\0'; sscanf(string,"%d",&pweight); comp=line[8]; string[0]=line[17]; string[1]=line[18]; string[2]=line[19]; string[3]=line[20]; string[4]='\0'; sscanf(string,"%d",&yrp); string[0]=line[21]; string[1]=line[22]; string[2]='\0'; sscanf(string,"%d",&mthp); string[0]=line[23]; string[1]=line[24]; string[2]='\0'; sscanf(string,"%d",&dayp); string[0]=line[25]; string[1]=line[26]; string[2]='\0'; sscanf(string,"%d",&hrp); string[0]=line[27]; string[1]=line[28]; string[2]='\0'; sscanf(string,"%d",&minp); for(i=0;i<5;++i)string[i]=line[i+29]; string[5]='\0'; sscanf(string,"%f",&secp); for(i=0;i<5;++i)string[i]=line[i+41]; string[5]='\0'; sscanf(string,"%f",&secs); for(i=0;i<3;++i)scom[i]=line[46+i]; scom[3]='\0'; string[0]=line[49]; string[1]='\0'; sscanf(string,"%d",&sweight); source=line[108]; secp/=100; secs/=100; doyp=doy(yrp,mthp,dayp); } /* if station is WWVB, IRG1, IRG2 continue to next card */ if(strcmp(station,"WWVB")==0)continue; if(strcmp(station,"IRG1")==0)continue; if(strcmp(station,"IRG2")==0)continue; /* if rtp data and not wanted continue to next card */ if(!rtp && (source=='R' || source=='P' || source=='O' || source=='W'))continue; /* store P pick if there is a good one */ if(pcomm[1]=='P' || pcomm[1]=='p'){ /* store P pick */ /* if 4 weight and not wanted go to S pick */ if(pweight>=4 && !fours)goto spick; /* if horizontal and not wanted go to S pick */ if((comp=='E' || comp=='N' || comp=='e' || comp=='n' || comp=='1' || comp=='2') && !phoriz)goto spick; seca=adjtime(yre,doye,hre,mine,sece,yrp,doyp,hrp,minp,secp); /* see if we have this station */ whichsta=-1; for(i=0;i=4 && !fours)continue; /* if vertical and not wanted go next phase card */ if(comp!='E' && comp!='N' && comp!='e' && comp!='n' && comp!='1' && comp!='2' && !svert)continue; /* checked it all so output s pick after adjusting time */ seca=adjtime(yre,doye,hre,mine,sece,yrp,doyp,hrp,minp,secs); /* see if we have this station */ whichsta=-1; for(i=0;i %s ; /bin/rm -f %s %s", nameout2,nameout3,nameout2,nameout3,nameout); system(line); } float adjtime(yre,doye,hre,mine,sece,yrp,doyp,hrp,minp,secp) /* subroutine to adjust pick time to origin time */ int yre,doye,hre,mine,yrp,doyp,hrp,minp; float sece,secp; { float seca; /* simplest case */ if(yre==yrp && doye==doyp && hre==hrp && mine==minp)seca=secp-sece; /* else if different minutes */ else if(yre==yrp && doye==doyp && hre==hrp){ seca=60.*(minp-mine)+secp-sece; } /* else if different hours */ else if(yre==yrp && doye==doyp){ minp+=60.*(hrp-hre); seca=60.*(minp-mine)+secp-sece; } /* else if different day of year */ else if(yre==yrp){ hrp+=24*(doyp-doye); minp+=60.*(hrp-hre); seca=60.*(minp-mine)+secp-sece; } /* the unlikely case of different years, but only 1 year difference */ /* I mean lets be real we are talking about first arrivals! */ else { doyp+=doy(yre,12,31); hrp+=24*(doyp-doye); minp+=60.*(hrp-hre); seca=60.*(minp-mine)+secp-sece; } return(seca); } #define LEN 100 #include getstr(query,str,insist) char query[],str[]; int insist; { int j,change; char line[LEN],line2[LEN]; start: ; change=0; j=0; while(!j){ printf("%s [%s]: ",query,str); if(fgets(line,LEN,stdin)==NULL)exit(0); if(strlen(line)>1){ j=sscanf(line,"%s",line2); if(j){ strcpy(str,line2); change=1; } } else j=1; } if(insist && !change)goto start; return(change); } getline(query,str,length,insist) char query[],str[]; int insist,length; { int j,change; int i; char line[LEN],line2[LEN]; start: ; change=0; j=0; while(!j){ printf("%s [%s]: ",query,str); if(fgets(line,LEN,stdin)==NULL)exit(0); j=strlen(line); if(j>1){ for(i=0;i1){ j=sscanf(line,"%f",&z); if(j){ *pfl=z; change=1; } } else j=1; } if(insist && !change)goto start; return(change); } getint(query,pint,insist) char query[]; int insist,*pint; { int j,change,i; char line[LEN]; start: ; change=0; j=0; while(!j){ printf("%s [%d]: ",query,*pint); if(fgets(line,LEN,stdin)==NULL)exit(0); if(strlen(line)>1){ j=sscanf(line,"%d",&i); if(j){ *pint=i; change=1; } } else j=1; } if(insist && !change)goto start; return(change); } getdouble(query,pdb,insist) char query[]; double *pdb; int insist; { int j,change; char line[LEN]; double z; start: ; change=0; j=0; while(!j){ printf("%s [%g]: ",query,*pdb); if(fgets(line,LEN,stdin)==NULL)exit(0); if(strlen(line)>1){ j=sscanf(line,"%lf",&z); if(j){ *pdb=z; change=1; } } else j=1; } if(insist && !change)goto start; return(change); } doy(year,mth,date) /* returns julian day */ /* if illegal mth or day returns 0 */ int year,mth,date; { static short ndays[12]={ 0,31,59,90,120,151,181,212,243,273,304,334 }; short days=0; short leap=0; if(mth<1 || mth >12) return (days); if(date<1 || date>31)return (days); days=ndays[mth-1]+date; if((year%4) == 0)leap=1; if((year%100) == 0)leap=0; if((year%1000) == 0)leap=1; if(leap && mth>2)days++; return (days); } mthchar(mth,month) /* returns 3 letter abbreviation of the month in month, leaves month unchanged if mth illegal */ /* also returns a 0 if succesful, -1 if mth illegal */ short mth; char month[]; { short i; short err= -1; static char mnths[37]= { "janfebmaraprmayjunjulaugsepoctnovdec" }; if(mth<1 || mth>12)return (err); err=0; i=3*(mth-1); month[0]=mnths[i]; month[1]=mnths[i+1]; month[2]=mnths[i+2]; return (err); }