Titan2Reader: sac.c Source File

00001 /*
00002  * This file is part of Titan2Reader
00003  *
00004  * Copyright (C) 2006, Sebastien Judenherc <sebastien.judenherc@agecodagis.com>
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307
00019  * USA
00020  *
00021  */
00022 
00023 #include <stdio.h>
00024 #include <ctype.h>
00025 #include <stdlib.h>
00026 #include <unistd.h>
00027 #include <time.h>
00028 #include <string.h>
00029 #include <assert.h>
00030 #include <sys/types.h>
00031 
00032 #include "common.h" 
00033 
00034 struct sac {                    
00035 float    delta,     depmin,    depmax,    scale,     odelta;    
00036 float    b,         e,         o,         a,         internal1; 
00037 float    t0,        t1,        t2,        t3,        t4;        
00038 float    t5,        t6,        t7,        t8,        t9;        
00039 float    f,         resp0,     resp1,     resp2,     resp3;     
00040 float    resp4,     resp5,     resp6,     resp7,     resp8;     
00041 float    resp9,     stla,      stlo,      stel,      stdp;      
00042 float    evla,      evlo,      evel,      evdp,      unused1;   
00043 float    user0,     user1,     user2,     user3,     user4;     
00044 float    user5,     user6,     user7,     user8,     user9;     
00045 float    dist,      az,        baz,       gcarc,     internal2; 
00046 float    internal3, depmen,    cmpaz,     cmpinc,    unused2;   
00047 float    unused3,   unused4,   unused5,   unused6,   unused7;   
00048 float    unused8,   unused9,   unused10,  unused11,  unused12;  
00049 int32_t    nzyear,    nzjday,    nzhour,    nzmin,     nzsec;     
00050 int32_t    nzmsec,    nvhdr,     internal5, internal6, npts;      
00051 int32_t    internal7, internal8, unused13,  unused14,  unused15;  
00052 int32_t    iftype,    idep,      iztype,    unused16,  iinst;     
00053 int32_t    istreg,    ievreg,    ievtyp,    iqual,     isynth;    
00054 int32_t    unused17,  unused18,  unused19,  unused20,  unused21;  
00055 int32_t    unused22,  unused23,  unused24,  unused25,  unused26;  
00056 int32_t    leven,     lpspol,    lovrok,    lcalda,    unused27;  
00057 char    kstnm[8],  kevnm[16];           
00058 char    khole[8],  ko[8],     ka[8];               
00059 char    kt0[8],    kt1[8],    kt2[8];              
00060 char    kt3[8],    kt4[8],    kt5[8];              
00061 char    kt6[8],    kt7[8],    kt8[8];              
00062 char    kt9[8],    kf[8],     kuser0[8];           
00063 char    kuser1[8], kuser2[8], kcmpnm[8];           
00064 char    knetwk[8], kdatrd[8], kinst[8];            
00065 };
00066 
00067 struct sac sacHeader[MAX_CHANNEL];
00068 FILE *sacChannelFile[MAX_CHANNEL] = {
00069         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00070         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00071         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00072         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00073         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00074 };
00075 char sacChannelFileName[MAX_CHANNEL][128];
00076 
00077 static struct sac sac_null = {
00078 -12345., -12345., -12345., -12345., -12345.,
00079 -12345., -12345., -12345., -12345., -12345.,
00080 -12345., -12345., -12345., -12345., -12345.,
00081 -12345., -12345., -12345., -12345., -12345.,
00082 -12345., -12345., -12345., -12345., -12345.,
00083 -12345., -12345., -12345., -12345., -12345.,
00084 -12345., -12345., -12345., -12345., -12345.,
00085 -12345., -12345., -12345., -12345., -12345.,
00086 -12345., -12345., -12345., -12345., -12345.,
00087 -12345., -12345., -12345., -12345., -12345.,
00088 -12345., -12345., -12345., -12345., -12345.,
00089 -12345., -12345., -12345., -12345., -12345.,
00090 -12345., -12345., -12345., -12345., -12345.,
00091 -12345., -12345., -12345., -12345., -12345.,
00092 -12345, -12345, -12345, -12345, -12345,
00093 -12345, -12345, -12345, -12345, -12345,
00094 -12345, -12345, -12345, -12345, -12345,
00095 -12345, -12345, -12345, -12345, -12345,
00096 -12345, -12345, -12345, -12345, -12345,
00097 -12345, -12345, -12345, -12345, -12345,
00098 -12345, -12345, -12345, -12345, -12345,
00099 -12345, -12345, -12345, -12345, -12345,
00100 { '-','1','2','3','4','5',' ',' ' },
00101 { '-','1','2','3','4','5',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ' },
00102 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00103 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00104 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00105 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00106 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00107 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00108 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00109 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00110 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00111 { '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
00112 { '-','1','2','3','4','5',' ',' ' }
00113 };
00114 
00115 
00116 #define IREAL   0 
00117 #define ITIME   1 
00118 #define IRLIM   2 
00119 #define IAMPH   3 
00120 #define IXY     4 
00121 #define IUNKN   5 
00122 #define IDISP   6 
00123 #define IVEL    7 
00124 #define IACC    8 
00125 #define IB      9 
00126 #define IDAY   10 
00127 #define IO     11 
00128 #define IA     12 
00129 #define IT0    13 
00130 #define IT1    14 
00131 #define IT2    15 
00132 #define IT3    16 
00133 #define IT4    17 
00134 #define IT5    18 
00135 #define IT6    19 
00136 #define IT7    20 
00137 #define IT8    21 
00138 #define IT9    22 
00139 #define IRADNV 23 
00140 #define ITANNV 24 
00141 #define IRADEV 25 
00142 #define ITANEV 26 
00143 #define INORTH 27 
00144 #define IEAST  28 
00145 #define IHORZA 29 
00146 #define IDOWN  30 
00147 #define IUP    31 
00148 #define ILLLBB 32 
00149 #define IWWSN1 33 
00150 #define IWWSN2 34 
00151 #define IHGLP  35 
00152 #define ISRO   36 
00153 #define INUCL  37 
00154 #define IPREN  38 
00155 #define IPOSTN 39 
00156 #define IQUAKE 40 
00157 #define IPREQ  41 
00158 #define IPOSTQ 42 
00159 #define ICHEM  43 
00160 #define IOTHER 44 
00161 #define IGOOD  45 
00162 #define IGLCH  46 
00163 #define IDROP  47 
00164 #define ILOWSN 48 
00165 #define IRLDTA 49 
00166 #define IVOLTS 50 
00167 #define INIV51 51 
00168 #define INIV52 52 
00169 #define INIV53 53 
00170 #define INIV54 54 
00171 #define INIV55 55 
00172 #define INIV56 56 
00173 #define INIV57 57 
00174 #define INIV58 58 
00175 #define INIV59 59 
00176 #define INIV60 60 
00177 
00178 double sac2unix(struct sac *sac0)
00179 {
00180         time_t t0;
00181         struct tm tm0;
00182         if (sac0->iztype==IUNKN) {
00183                 PrintError("iztype is IUNKN\n");
00184                 return(-1.0);
00185         }
00186         tm0.tm_mday=sac0->nzjday;
00187         tm0.tm_year=sac0->nzyear-1900;
00188         tm0.tm_hour=sac0->nzhour;
00189         tm0.tm_min=sac0->nzmin;
00190         tm0.tm_sec=sac0->nzsec;
00191         tm0.tm_mon=0;
00192         tm0.tm_isdst=0;
00193         t0=mktime(&tm0);
00194         if (t0<0) {
00195                 PrintError("mktime error %d %d %d %d %d\n",
00196                                 tm0.tm_mday,
00197                                 tm0.tm_year,
00198                                 tm0.tm_hour,
00199                                 tm0.tm_min,
00200                                 tm0.tm_sec,
00201                                 tm0.tm_mon,
00202                                 tm0.tm_isdst);
00203                 return((double)t0);
00204         }
00205         return((double)t0+0.001*(double)sac0->nzmsec);
00206 }
00207 
00208 void unix2sac(struct sac *sac0,double temps)
00209 {
00210         time_t t0=(time_t)temps;
00211         struct tm tm0;
00212         if (t0<0) PrintError("mktime error\n");
00213         memcpy(&tm0,gmtime(&t0),sizeof(struct tm));
00214         sac0->nzjday=tm0.tm_yday+1;
00215         sac0->nzyear=tm0.tm_year+1900;
00216         sac0->nzhour=tm0.tm_hour;
00217         sac0->nzmin=tm0.tm_min;
00218         sac0->nzsec=tm0.tm_sec;
00219         sac0->nzmsec=(int32_t)(1000.*(temps-(double)t0));
00220 }
00221 
00222 int32_t defsack(char * field)
00223 {
00224         if (!field) {
00225                 return(0); /*...*/
00226         }
00227         return(memcmp(field,"-12345",6));
00228 }
00229 
00233 FILE* OpenSacChannel(int channel)
00234 {
00235         if (sacChannelFile[channel]==NULL)
00236         {
00237                 int fd;
00238                 strcpy(sacChannelFileName[channel],"titanfileXXXXXX");
00239                 fd=OpenTempFile(sacChannelFileName[channel]);
00240                 assert((sacChannelFile[channel]=fdopen(fd,"wb"))!=NULL);
00241                 setvbuf(sacChannelFile[channel],NULL,_IOFBF,FILE_BUFFER_SIZE);
00242 
00243                 if (sacChannelFile[channel])
00244                 {
00245                         // initialized a new generic header
00246                         InitInfoHeader(channel);
00247                         // initiliaze a new sac header
00248                         memcpy(&(sacHeader[channel]),&sac_null,sizeof(struct sac));
00249                         sacHeader[channel].lovrok=1;
00250                         sacHeader[channel].lcalda=1;
00251                         sacHeader[channel].iftype=ITIME;
00252                         sacHeader[channel].leven=1;
00253                         sacHeader[channel].npts=0;
00254                         sacHeader[channel].b=0.;
00255                         sacHeader[channel].nvhdr=6;
00256                         sacHeader[channel].iztype=IB;
00257                         sacHeader[channel].delta=(double)HWConfig.srExp[channel].div/
00258                                 (double)HWConfig.srExp[channel].base;
00259                         // and write the header to the file
00260                         fwrite(&sacHeader[channel],sizeof(struct sac),1,sacChannelFile[channel]);
00261                 }
00262                 else
00263                 {
00264                         PrintError("unable to open %s: %m\n",sacChannelFileName[channel]);
00265                 }
00266         }
00267         return(sacChannelFile[channel]);
00268 }
00269 
00273 void CloseSacFile(int channel)
00274 {
00275         char line[128];
00276         if (sacChannelFile[channel]==NULL)
00277         {
00278                 return;
00279         }
00280         // update the internal information header and the sac start time
00281         unix2sac(&(sacHeader[channel]),UpdateInfoHeader(channel));
00282         // then copy the appropriate fields to the SAC header
00283         // kstnm
00284         snprintf(line,8,"%-8s",infoHeader[channel].station);
00285         memcpy(sacHeader[channel].kstnm,line,8);
00286         // kcmpnm
00287         if (strlen(infoHeader[channel].channel)>0)
00288         {
00289                 snprintf(line,8,"%-8s",infoHeader[channel].channel);
00290         }
00291         else
00292         {
00293                 snprintf(line,8,"%-8d",channel);
00294         }
00295         memcpy(sacHeader[channel].kcmpnm,line,8);
00296 
00297         // update station information
00298         sacHeader[channel].stla=oldHWConfig.avgPos.lat;
00299         sacHeader[channel].stlo=oldHWConfig.avgPos.lon;
00300         sacHeader[channel].stel=oldHWConfig.avgPos.elev;
00301 
00302         // update source information
00303         sacHeader[channel].evla=event.lat;
00304         sacHeader[channel].evlo=event.lon;
00305         sacHeader[channel].evdp=event.dep;
00306         sacHeader[channel].o=event.OT-infoHeader[channel].correctedTime;
00307 
00308         // rewrite the header
00309         fseek(sacChannelFile[channel],0,SEEK_SET);
00310         fwrite(&sacHeader[channel],sizeof(struct sac),1,sacChannelFile[channel]);
00311         fclose(sacChannelFile[channel]);
00312         sacChannelFile[channel]=NULL;
00313         strcpy(line,FileName("", ".SAC", channel, fileNameOptions));
00314         if (noWrite==0)
00315         {
00316                 Rename(sacChannelFileName[channel],line);
00317         }
00318         PrintLog("OUTPUT: %s %d samples (%.2fs)\n",line,
00319                         infoHeader[channel].npts,
00320                         (float)(infoHeader[channel].npts*
00321                                 infoHeader[channel].sr.div)/
00322                                         (float)infoHeader[channel].sr.base);
00323         unlink(sacChannelFileName[channel]);
00324 }
00325 
00328 void AbortSacFile(int channel)
00329 {
00330         PrintDebug("abort for %s\n",sacChannelFileName[channel]);
00331         if (sacChannelFile[channel])
00332         {
00333                 fclose(sacChannelFile[channel]);
00334                 sacChannelFile[channel]=NULL;
00335                 unlink(sacChannelFileName[channel]);
00336         }
00337 }
00338 
00341 void AddSacData(int channel)
00342 {
00343         float fbuff[2000];
00344         int i;
00345         if ((HWConfig.srObs[channel].npts==0)||(HWConfig.srExp[channel].npts==0))
00346         {
00347                 return;
00348         }
00349         if (OpenSacChannel(channel))
00350         {
00351                 for (i=0; i<HWConfig.srObs[channel].npts; i++)
00352                 {
00353                         fbuff[i]=(float)HWConfig.data[channel][i];
00354                 }
00355                 if (privateDebug)
00356                 {
00357                         if (channel==2)
00358                         {
00359                                 fbuff[0]=9999.;
00360                         }
00361                 }
00362                 if (noWrite==0)
00363                 {
00364                         if (fwrite(fbuff,4,  
00365                                                 HWConfig.srObs[channel].npts,sacChannelFile[channel])
00366                                         !=HWConfig.srObs[channel].npts)
00367                         {
00368                                 PrintError("unable to write to %s: %m\n",sacChannelFileName[channel]);
00369                                 AbortSacFile(channel);
00370                         }
00371                 }
00372                 infoHeader[channel].npts+=HWConfig.srObs[channel].npts;
00373                 sacHeader[channel].npts+=HWConfig.srObs[channel].npts;
00374         }
00375         if (test0==1)
00376         {
00377                 int nz=0;
00378                 for (i=0; i<HWConfig.srObs[channel].npts; i++)
00379                 {
00380                         if (HWConfig.data[channel][i]==0)
00381                         {
00382                                 nz++;
00383                         }
00384                         else
00385                         {
00386                                 nz=0;
00387                         }
00388                 }
00389                 if (nz>5)
00390                 {
00391                         PrintError("%d consecutive null data found @ %s [%ld:%s]\n",
00392                                         StrTime(&currentTime),
00393                                         currentFile,
00394                                         SafeFtell(titan2File));
00395                         PrintLog("%d consecutive null data found @ %s [%ld:%s]\n",
00396                                         StrTime(&currentTime),
00397                                         currentFile,
00398                                         SafeFtell(titan2File));
00399 
00400                 }
00401         }
00402 }
00403 
00409 void SetSource(float lat, float lon, float dep, double OT)
00410 {
00411         event.lat=lat;
00412         event.lon=lon;
00413         event.dep=dep;
00414         event.OT=OT;
00415 }