Titan2Reader: seed.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 #include <netinet/in.h>
00032 
00033 #include "common.h" 
00034 
00035 FILE *mseedChannelFile[MAX_CHANNEL] = {
00036         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00037         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00038         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00039         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00040         NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,
00041 };
00042 char mseedChannelFileName[MAX_CHANNEL][128];
00043 
00044 // data encoding, default is 3 (32bits integers)
00045 int seedDataFormat=3;
00046 
00047 double staticDt[MAX_CHANNEL];
00048 
00049 #define US2S 0.23283064365386962891e-9
00050 
00051 #define FSDH_SIZE 48
00052 
00053 #define BLOCKLEN (1<<(9+seedBlockSizePowerOf2))
00054 
00055 #define NETBO
00056 
00057 #ifdef NETBO
00058 #       define ToNBOs(x) htons(x)
00059 #       define ToNBOl(x) htonl(x)
00060 #       define ToHBOs(x) ntohs(x)
00061 #       define ToHBOl(x) ntohl(x)
00062 #else
00063 #       define ToNBOs(x) (x)
00064 #       define ToNBOl(x) (x)
00065 #       define ToHBOs(x) (x)
00066 #       define ToHBOl(x) (x)
00067 #endif
00068 
00069 
00070 /*  SEED FSDH: FixedSeedDataHeader_b1000_b1001 Section Data Header */
00071 
00072 #pragma pack(1)
00073 struct SeedTime
00074 {
00075         int16_t        year;
00076         int16_t        day;
00077         char            hour;
00078         char            min;
00079         char            sec;
00080         char            unused;
00081         int16_t        frac10000;
00082 }; //__attribute__((packed));
00083 
00084 #pragma pack(1)
00085 struct FixedSeedDataHeader_b1000_b1001
00086 {
00087         char            sequence[6];
00088         char            dhquality;
00089         char            reserved;
00090         char            station[5];
00091         char            location[2];
00092         char            channel[3];
00093         char            network[2];
00094         struct SeedTime beginTime;
00095         int16_t        numSamples;
00096         int16_t         sampleRateFactor;
00097         int16_t         sampleRateMultiplier;
00098         char            activityFlags;
00099         char            ioFlags;
00100         char            dataQualityFlags;
00101         char            numBlockettes;
00102         int32_t         timeCorrection;
00103         int16_t        beginData;
00104         int16_t        firstBlockette;
00105         // blockette 1000
00106         int16_t b1000_type;
00107         int16_t b1000_next; // 0 si derni?re
00108         char b1000_format;
00109         char b1000_order;
00110         char b1000_length;
00111         char b1000_reserved;
00112         // blockette 1001
00113         int16_t b1001_type;
00114         int16_t b1001_next;
00115         unsigned char b1001_timingQ;
00116         unsigned char b1001_microSecond;
00117         unsigned char b1001_reserved;
00118         unsigned char b1001_frameCount;
00119 }; // __attribute__((packed));
00120 
00121 #define HEADEROFFSET 16
00122 
00123 int ParseSeedOptions(char *cmd)
00124 {
00125         char *s=cmd;
00126         int next,error;
00127         if (*s==0)
00128         {
00129                 return(0);
00130         }
00131         if (*s!=',')
00132         {
00133                 return(-1);
00134         }
00135         while (*s)
00136         {
00137                 s++;
00138                 next=0;
00139                 error=1;
00140                 if (strncmp(s,"f=",2)==0)
00141                 {
00142                         if (sscanf(s+2,"%d%n",&seedDataFormat,&next)!=1)
00143                         {
00144                                 break;
00145                         }
00146                         // known: 3, 10
00147                         if ((seedDataFormat!=3)&&(seedDataFormat!=10))
00148                         {
00149                                 PrintError("SEED format %d not supported\n",seedDataFormat);
00150                                 return(-1);
00151                         }
00152                         PrintDebug("format=%d\n",seedDataFormat);
00153                         error=0;
00154                         next+=2;
00155                 }
00156                 if (strncmp(s,"b=",2)==0)
00157                 {
00158                         if (sscanf(s+2,"%d%n",&seedBlockSizePowerOf2,&next)!=1)
00159                         {
00160                                 break;
00161                         }
00162                         // 0:512, 1:1K 2:2K 3:4K
00163                         if ((seedBlockSizePowerOf2<0)||(seedBlockSizePowerOf2>3))
00164                         {
00165                                 PrintError("SEED block length out of range\n");
00166                                 return(-1);
00167                         }
00168                         PrintDebug("blocklength=%d bytes\n",BLOCKLEN);
00169                         error=0;
00170                         next+=2;
00171                 }
00172                 if (error)
00173                 {
00174                         PrintError("could not parse mseed option string\n");
00175                         return(-1);
00176                 }
00177                 s+=next;
00178         }
00179         return(0);
00180 }
00181 
00182 struct SeedTime *TimeToSeedTime(double T0, struct SeedTime *sd)
00183 {
00184         time_t tt;
00185         struct tm tm;
00186         int frac;
00187         tt=(time_t)T0;
00188         memcpy(&tm,gmtime(&tt),sizeof(struct tm));
00189         sd->year=ToNBOs(tm.tm_year+1900);
00190         sd->day=ToNBOs(tm.tm_yday+1);
00191         sd->hour=tm.tm_hour;
00192         sd->min=tm.tm_min;
00193         sd->sec=tm.tm_sec;
00194         sd->unused=0;
00195         frac=(int)(10000.*(T0-(double)tt));
00196         if (frac>9999)
00197         {
00198                 sd->sec++;
00199                 frac-=10000;
00200         }
00201         sd->frac10000=ToNBOs((int16_t)(frac));
00202         return(sd);
00203 }
00204 
00205 double SeedTimeToTime(struct SeedTime *sd)
00206 {
00207         struct tm tm;
00208         memset(&tm,0,sizeof(struct tm));
00209         tm.tm_year=(int)ToHBOs(sd->year)-1900;
00210         tm.tm_mon=0;
00211         tm.tm_mday=(int)ToHBOs(sd->day);
00212         tm.tm_hour=sd->hour;
00213         tm.tm_min=sd->min;
00214         tm.tm_sec=sd->sec;
00215         return((double)mktime(&tm)+0.0001*(double)ToHBOs(sd->frac10000));
00216 }
00217 
00218 
00219 void UpdateFSDH(struct FixedSeedDataHeader_b1000_b1001 *fsdh,
00220                 int seq, int ndata, double CT, double correction,
00221                 char *station, char *location,
00222                 char *channel, char *network)
00223 {
00224         char tmp[10];
00225 
00226         // voir "Seed Manual v2.4", chap 8
00227 
00228         sprintf(tmp, "%06d", seq);
00229         memcpy(fsdh->sequence, tmp, 6);
00230 
00231         // blk 1000
00232         fsdh->b1000_type=ToNBOs(1000);
00233         fsdh->b1000_next=ToNBOs(48+8); // 48=FSDH, 8=1000
00234         fsdh->b1000_format=seedDataFormat; // 2=int 24 bits, 3=int 32 bits, 10=steim1
00235 #ifdef NETBO
00236         fsdh->b1000_order=1;  // 1=Big Endian
00237 #else
00238         fsdh->b1000_order=0;  // 0=Little Endian
00239 #endif
00240         fsdh->b1000_length=9+seedBlockSizePowerOf2;
00241 
00242         // blk 1001
00243         fsdh->b1001_type=ToNBOs(1001);
00244         fsdh->b1001_next=0; // last one
00245         fsdh->b1001_timingQ=0;
00246         // fsdh->b1001_microSecond later
00247         fsdh->b1001_reserved=32;
00248 
00249         // FSDH
00250         fsdh->dhquality='D';
00251         fsdh->reserved=32;
00252 
00253         // station
00254         sprintf(tmp,"%-5s",station);
00255         memcpy(fsdh->station,tmp,5);
00256         // location
00257         sprintf(tmp,"%-2s",location);
00258         memcpy(fsdh->location,tmp,2);
00259         // location
00260         sprintf(tmp,"%-3s",channel);
00261         memcpy(fsdh->channel,tmp,3);
00262         // network
00263         sprintf(tmp,"%-2s",network);
00264         memcpy(fsdh->network,tmp,2);
00265 
00266         // beginTime
00267         if ((CT<1.e10)&&(CT>-1.e10))
00268         {
00269                 TimeToSeedTime(CT,&(fsdh->beginTime));
00270                 fsdh->b1001_microSecond=(char)(100.*((CT*1.e4)-(double)((int64_t)(CT*1.e4))));
00271         }
00272 
00273         // nb d'échantillons
00274         if (ndata>0)
00275         {
00276                 fsdh->numSamples=ToNBOs((int16_t)ndata);
00277         }
00278 
00279         // activity flags: bit1=1(correction appliqu?e)
00280         fsdh->activityFlags=1<<1;
00281 
00282         // io and clock flags
00283         fsdh->ioFlags=0; // bit5=clock locked
00284 
00285         // data quality flags
00286         fsdh->dataQualityFlags=(1<<7); // bit7=possible problème de temps, corrigé plus tard
00287 
00288         // nb de blockettes ? suivre
00289         fsdh->numBlockettes=2; // blockettes 1000 and 1001
00290 
00291         // la correction de temps
00292         if ((correction<1.e10)&&(correction>-1.e10))
00293         {
00294                 // mis à jour en deuxième passe
00295                 fsdh->timeCorrection=ToNBOl((int32_t)(correction*1.e4));
00296                 fsdh->b1001_timingQ=100; // now we know the time
00297                 fsdh->dataQualityFlags=0;
00298         }
00299 
00300         // d?but des donn?es
00301         if (seedDataFormat==10)
00302         {
00303                 fsdh->beginData=ToNBOs(HEADEROFFSET*4);
00304                 fsdh->b1001_frameCount=(unsigned char)(BLOCKLEN-4*HEADEROFFSET)/64;
00305         }
00306         else
00307         {
00308                 fsdh->beginData=ToNBOs(sizeof(struct FixedSeedDataHeader_b1000_b1001));
00309                 fsdh->b1001_frameCount=(unsigned char)(BLOCKLEN-sizeof(struct FixedSeedDataHeader_b1000_b1001))/64;
00310         }
00311 
00312         // first blockette
00313         fsdh->firstBlockette=ToNBOs(FSDH_SIZE);
00314 }
00315 
00316 
00320 FILE* OpenMSeedChannel(int channel)
00321 {
00322         // ask here for the namùing options
00323         fileNameOptions|=(ADD_LOCATION_CODE|ADD_NETWORK_CODE);
00324         if (mseedChannelFile[channel]==NULL)
00325         {
00326                 int fd;
00327                 strcpy(mseedChannelFileName[channel],"titanfileXXXXXX");
00328                 fd=OpenTempFile(mseedChannelFileName[channel]);
00329                 assert((mseedChannelFile[channel]=fdopen(fd,"w+b"))!=NULL);
00330                 setvbuf(mseedChannelFile[channel],NULL,_IOFBF,FILE_BUFFER_SIZE);
00331 
00332                 if (mseedChannelFile[channel])
00333                 {
00334                         // initialized a new generic header
00335                         InitInfoHeader(channel);
00336                 }
00337                 else
00338                 {
00339                         PrintError("unable to open %s: %m\n",mseedChannelFileName[channel]);
00340                 }
00341         }
00342         return(mseedChannelFile[channel]);
00343 }
00344 
00348 void CloseMSeedInt32File(int channel)
00349 {
00350         char line[128];
00351         char block[BLOCKLEN];
00352         struct FixedSeedDataHeader_b1000_b1001 *fsdh;
00353         int modlen=0,tnd=0,nd=0;
00354         int seq=0;
00355         int rs,total,noTimeWarning=0;
00356         if (mseedChannelFile[channel]==NULL)
00357         {
00358                 return;
00359         }
00360         // update the header information
00361         UpdateInfoHeader(channel);
00362         infoHeader[channel].correctionMethod=INTERPOLATION_CORRECTION;
00363 
00364         modlen=(BLOCKLEN-sizeof(struct FixedSeedDataHeader_b1000_b1001))/4;
00365         // reread all the fsdh and update each of them (time and number of samples)
00366         fsdh=(struct FixedSeedDataHeader_b1000_b1001*)block;
00367         fseek(mseedChannelFile[channel],0,SEEK_END);
00368         total=ftell(mseedChannelFile[channel]);
00369         fseek(mseedChannelFile[channel],0,SEEK_SET);
00370         while ((rs=fread(block,1,BLOCKLEN,mseedChannelFile[channel]))>
00371                         sizeof(struct FixedSeedDataHeader_b1000_b1001))
00372         {
00373                 double T,dt;
00374                 int age;
00375                 T=SeedTimeToTime(&(fsdh->beginTime));
00376                 T+=1.e-6*(double)fsdh->b1001_microSecond;
00377                 T-=staticDt[channel];
00378                 if (GetInterpolatedShift(T, &dt, &age, infoHeader[channel].station, NULL)
00379                                 ==NO_TIME_CORRECTION)
00380                 {
00381                         dt=infoHeader[channel].otshift;
00382                         PrintWarning("no time interpolation for block %d, using observed shift: %.3fs\n",seq,dt);
00383                 }
00384                 T-=dt;
00385                 if (infoHeader[channel].correctionMethod==NO_TIME_CORRECTION)
00386                 {
00387                         dt=1.e12;
00388                         noTimeWarning=1;
00389                 }
00390                 PrintDebug("T=%s, dt=%f seq=%d age=%d\n",time2str(T),dt,seq,age);
00391                 if (tnd+modlen>=infoHeader[channel].npts)
00392                 {
00393                         nd=infoHeader[channel].npts-tnd;
00394                 }
00395                 else
00396                 {
00397                         nd=modlen;
00398                 }
00399 
00400                 UpdateFSDH(
00401                                 fsdh,
00402                                 seq,
00403                                 nd, // ndata: modlen
00404                                 T, // corected time
00405                                 dt, // time correction
00406                                 infoHeader[channel].station,
00407                                 infoHeader[channel].location,
00408                                 infoHeader[channel].channel,
00409                                 infoHeader[channel].network);
00410                 fseek(mseedChannelFile[channel],seq*BLOCKLEN,SEEK_SET);
00411                 fwrite(block,BLOCKLEN,1,mseedChannelFile[channel]);
00412                 fflush(mseedChannelFile[channel]);
00413                 rs=ftell(mseedChannelFile[channel]);
00414                 if (rs>=total) // why the fread is not able to see the end of file ???
00415                 {
00416                         break;
00417                 }
00418                 PrintDebug("pos=%d/%d\n",rs,total);
00419                 tnd+=nd;
00420                 seq++;
00421         }
00422 
00423         fclose(mseedChannelFile[channel]);
00424         mseedChannelFile[channel]=NULL;
00425         strcpy(line,FileName("", ".mseed", channel, fileNameOptions));
00426         if (noWrite==0)
00427         {
00428                 Rename(mseedChannelFileName[channel],line);
00429         }
00430         PrintLog("OUTPUT: %s %d samples (%.2fs)\n",line,infoHeader[channel].npts,
00431                         (float)infoHeader[channel].npts/(float)infoHeader[channel].sr.npts);
00432         if (noTimeWarning==1)
00433         {
00434                 PrintWarning("The time correction could not be interpolated "
00435                                 "for some blocks of this file. "
00436                                 "Make sure that the time "
00437                                 "correction file exists. "
00438                                 "A preliminary pass with -time may solve the "
00439                                 "problem.\n");
00440         }
00441         unlink(mseedChannelFileName[channel]);
00442 }
00443 
00446 void AbortMSeedFile(int channel)
00447 {
00448         PrintDebug("abort for %s\n",mseedChannelFileName[channel]);
00449         if (mseedChannelFile[channel])
00450         {
00451                 fclose(mseedChannelFile[channel]);
00452                 mseedChannelFile[channel]=NULL;
00453                 unlink(mseedChannelFileName[channel]);
00454         }
00455 }
00456 
00459 void AddMSeedDataInt32(int channel)
00460 {
00461         int i,seq;
00462         int32_t sample;
00463         struct FixedSeedDataHeader_b1000_b1001 fsdh;
00464         int modlen;
00465         seq=0;
00466         if ((currentTime.usecond>>32)<1)
00467         {
00468                 PrintWarning("incomplete time frame (%lld : %lld): end of input ?\n",
00469                                 currentTime.usecond,currentTime.usecond>>32);
00470                 return;
00471         }
00472         if ((HWConfig.srObs[channel].npts==0)||(HWConfig.srExp[channel].npts==0))
00473         {
00474                 return;
00475         }
00476         // n sample
00477         modlen=(BLOCKLEN-sizeof(struct FixedSeedDataHeader_b1000_b1001))/4;
00478         if (OpenMSeedChannel(channel))
00479         {
00480                 // for each sample
00481                 for (i=0; i<HWConfig.srObs[channel].npts; i++)
00482                 {
00483                         // if this sample is the first of a new block
00484                         if ((infoHeader[channel].npts%modlen)==0)
00485                         {
00486                                 // write an empty fsdh (time=uncorrected)
00487                                 double ut,dt;
00488                                 int ch=PrimaryChannel(channel);
00489                                 ut=US2S*(double)(currentTime.usecond-
00490                                                 ((int64_t)((HWConfig.srExp[ch].npts)*HWConfig.srExp[ch].div)<<32)/
00491                                                 (int64_t)HWConfig.srExp[channel].base);
00492                                 dt=(double)((i+1)*HWConfig.srExp[channel].div)/
00493                                         (double)HWConfig.srExp[channel].base;
00494                                 ut+=dt;
00495 
00496                                 memset(&fsdh,0,sizeof(struct FixedSeedDataHeader_b1000_b1001));
00497 
00498                                 UpdateFSDH(
00499                                                 &fsdh,
00500                                                 seq,
00501                                                 modlen, // ndata: modlen
00502                                                 ut, // uncorrected time
00503                                                 1.e12, // time correction
00504                                                 HWConfig.fieldIdent,
00505                                                 HWConfig.locationName[channel],
00506                                                 HWConfig.channelName[channel],
00507                                                 HWConfig.networkName);
00508                                 // sampling rate
00509                                 fsdh.sampleRateFactor=ToNBOs(HWConfig.srExp[channel].base);
00510                                 fsdh.sampleRateMultiplier=ToNBOs(-HWConfig.srExp[channel].div);
00511 
00512                                 fwrite(&(fsdh),sizeof(struct FixedSeedDataHeader_b1000_b1001),
00513                                                 1,mseedChannelFile[channel]);
00514 
00515                                 seq++;
00516                         }
00517                         // write the sample
00518                         sample=ToNBOl(HWConfig.data[channel][i]);
00519                         fwrite(&sample,4,1,mseedChannelFile[channel]);
00520                         // increase the number of samples
00521                         infoHeader[channel].npts++;
00522                 // end for
00523                 }
00524         }
00525 }
00526 
00529 void AddMSeedDataSteim1(int channel)
00530 {
00531         int i;
00532         double T0,T;
00533         int ch=PrimaryChannel(channel);
00534         if ((currentTime.usecond>>32)<1)
00535         {
00536                 PrintWarning("incomplete time frame (%lld : %lld): end of input ?\n",
00537                                 currentTime.usecond,currentTime.usecond>>32);
00538                 return;
00539         }
00540         if (
00541                         (HWConfig.srObs[channel].npts==0)||
00542                         (HWConfig.srObs[channel].npts!=HWConfig.srExp[channel].npts)||
00543                         (HWConfig.srExp[channel].npts==0))
00544         {
00545                 PrintWarning("skipping partial superframe\n");
00546                 return;
00547         }
00548         // compute the time of the first sample
00549         T0=US2S*(double)(currentTime.usecond-
00550                         ((int64_t)((HWConfig.srExp[ch].npts)*HWConfig.srExp[ch].div)<<32)/
00551                         (int64_t)HWConfig.srExp[channel].base);
00552         if (OpenMSeedChannel(channel))
00553         {
00554                 // for each sample
00555                 for (i=0; i<HWConfig.srObs[channel].npts; i++)
00556                 {
00557                         // correct from the number of samples already written
00558                         T=T0+(double)((i+1)*HWConfig.srExp[channel].div)/
00559                                 (double)HWConfig.srExp[channel].base;
00560                         // write the time and the sample
00561                         fwrite(&T,sizeof(double),1,mseedChannelFile[channel]);
00562                         fwrite(&(HWConfig.data[channel][i]),sizeof(int32_t),1,mseedChannelFile[channel]);
00563                         // increase the number of samples
00564                         infoHeader[channel].npts++;
00565                 // end for
00566                 }
00567         }
00568 }
00569 
00572 void CloseMSeedSteim1File(int channel)
00573 {
00574         int32_t position,lastPosition,previous,lastPrevious,first;
00575         int32_t nvalues[4]={ 4, 4, 2, 1 };
00576         int32_t d,w[4],k,s;
00577         int32_t wk=0;
00578         int32_t c;
00579         int32_t nsubframe,wndx,ndata;
00580         int32_t seq=0;
00581         double t;
00582         FILE *fout;
00583         struct FixedSeedDataHeader_b1000_b1001 fsdh;
00584         int32_t frame[4096]; // 4096 is the maximum length for SEED frames
00585         char outputFileName[128];
00586         int noTimeWarning=0;
00587         int32_t totalndata=0;
00588         if (mseedChannelFile[channel]==NULL)
00589         {
00590                 return;
00591         }
00592 
00593         // init
00594         lastPrevious=previous=0;
00595         lastPosition=position=0;
00596         nsubframe=0;
00597         first=0xffffffff;
00598         k=c=0;
00599         memset(frame,0,BLOCKLEN);
00600         // first word index is 3 (c00, X0, XN ...)
00601         wndx=2;
00602         ndata=0;
00603         // close the file and reopen it for reading
00604         fclose(mseedChannelFile[channel]);
00605         mseedChannelFile[channel]=fopen(mseedChannelFileName[channel],"rb");
00606         if (mseedChannelFile[channel]==NULL)
00607         {
00608                 PrintError("could not open temp file for reading, abort\n");
00609                 unlink(mseedChannelFileName[channel]);
00610                 return;
00611         }
00612         // open the output for writing
00613         UpdateInfoHeader(channel);
00614         strcpy(outputFileName,FileName("", ".mseed", channel, fileNameOptions));
00615         fout=fopen(outputFileName,"wb");
00616         if (fout==NULL)
00617         {
00618                 PrintError("could not open %s for writing, abort\n",outputFileName);
00619                 fclose(mseedChannelFile[channel]);
00620                 unlink(mseedChannelFileName[channel]);
00621                 return;
00622         }
00623 
00624         // do
00625         do {
00626                 // seek to the current position
00627                 if (0!=fseek(mseedChannelFile[channel],position,SEEK_SET))
00628                 {
00629                         PrintError("%m\n");
00630                 }
00631                 // read the time and the sample
00632                 if (1!=fread(&t,sizeof(double), 1,mseedChannelFile[channel]))
00633                 {
00634                         PrintError("%m\n");
00635                 }
00636                 if (1!=fread(&s,sizeof(int32_t),1,mseedChannelFile[channel]))
00637                 {
00638                         PrintError("%m\n");
00639                 }
00640                 // update the first value if necessary
00641                 if (first==0xffffffff)
00642                 {
00643                         first=s;
00644                 }
00645                 // position=next position
00646                 position+=(sizeof(double)+sizeof(int32_t));
00647                 // compute delta
00648                 d=s-previous;
00649                 // update previous
00650                 previous=s;
00651                 // update w[k]
00652                 w[k]=d;
00653                 // increase k
00654                 k++;
00655                 // update c
00656                 switch (c)
00657                 {
00658                         case 0:
00659                         case 1:
00660                                 if ((d>=-128)&&(d<=127))
00661                                 {
00662                                         c=1;
00663                                         break;
00664                                 }
00665                         case 2:
00666                                 if ((d>=-32768)&&(d<=32767))
00667                                 {
00668                                         c=2;
00669                                         break;
00670                                 }
00671                                 c=3;
00672                         case 3:
00673                                 break;
00674                         default:
00675                                 PrintError("should never be there!!!!\n");
00676                                 abort();
00677                 }
00678                 // if we are about to run out of data
00679                 if (infoHeader[channel].npts-totalndata<nvalues[c])
00680                 {
00681                         // skip compression
00682                         c=3; // short cut
00683                 // endif
00684                 }
00685                 // if k>nvalues[c]
00686                 if (k>nvalues[c])
00687                 {
00688                         // leave c
00689                         // roll back the previous samples
00690                         previous=lastPrevious;
00691                         // reset k
00692                         k=0;
00693                         // reset position to lastPosition
00694                         position=lastPosition;
00695                         // skip the end of the loop
00696                         continue;
00697                         // endif
00698                 }
00699                 // if k==nvalues[c]
00700                 if (k>=nvalues[c])
00701                 {
00702                         assert(k==nvalues[c]);
00703                         // lastPosition=position
00704                         lastPosition=position;
00705                         // lastprevious=previous
00706                         lastPrevious=previous;
00707                         // compute the word
00708                         switch (c)
00709                         {
00710                                 case 1:
00711                                         wk =(w[0]&0xff)<<24;
00712                                         wk|=(w[1]&0xff)<<16;
00713                                         wk|=(w[2]&0xff)<<8;
00714                                         wk|=(w[3]&0xff);
00715                                         break;
00716                                 case 2:
00717                                         wk =(w[0]&0xffffL)<<16;
00718                                         wk|=(w[1]&0xffffL);
00719                                         break;
00720                                 case 3:
00721                                         wk=w[0];
00722                                         break;
00723                                 default:
00724                                         PrintError("unexepected c=%d\n",c);
00725                         }
00726                         // update the subframe
00727                         /*
00728                          PrintError("c=%d k=%d wk=%#0x ndx=%d nw=%d nsubframe=%d\n",c,k,wk,
00729                                         HEADEROFFSET+16*nsubframe+1+wndx,
00730                                         wndx,nsubframe);
00731                                         */
00732                         frame[HEADEROFFSET+16*nsubframe]|=(c<<(2*(14-wndx)));
00733                         frame[HEADEROFFSET+16*nsubframe+1+wndx]=wk;
00734                         // increase the number of words
00735                         wndx++;
00736                         // increase the number of sample of the subframe
00737                         ndata+=k;
00738                         totalndata+=k;
00739                         // reset c,k,w
00740                         c=k=0;
00741                         // if the subframe is full
00742                         if (wndx>=15)
00743                         {
00744                                 // increase the number of subframes
00745                                 nsubframe++;
00746                                 // reset the number of words
00747                                 wndx=0;
00748                         // endif
00749                         }
00750                 // endif
00751                 }
00752                 // if the frame is full or end of data reached
00753                 if ((4*(HEADEROFFSET+16*nsubframe)>=BLOCKLEN)||(totalndata==infoHeader[channel].npts))
00754                 {
00755                         double dt;
00756                         int32_t age;
00757                         int i;
00758                         // compute the time of the first sample of the block
00759                         t-=(double)((ndata-1)*HWConfig.srExp[channel].div)/
00760                                 (double)HWConfig.srExp[channel].base;
00761                         t-=staticDt[channel];
00762                         // compute the correction
00763                         if (GetInterpolatedShift(t, &dt, &age, infoHeader[channel].station, NULL)
00764                                         ==NO_TIME_CORRECTION)
00765                         {
00766                                 dt=infoHeader[channel].otshift;
00767                                 PrintWarning("no time interpolation for block %d, using observed shift: %.3fs\n",seq,dt);
00768                         }
00769                         t-=dt;
00770                         if (infoHeader[channel].correctionMethod==NO_TIME_CORRECTION)
00771                         {
00772                                 dt=1.e12;
00773                                 noTimeWarning=1;
00774                         }
00775                         // create the correct FSDH
00776                         UpdateFSDH(
00777                                         &fsdh,
00778                                         seq,
00779                                         ndata, // ndata
00780                                         t, // corected time
00781                                         dt, // time correction
00782                                         infoHeader[channel].station,
00783                                         infoHeader[channel].location,
00784                                         infoHeader[channel].channel,
00785                                         infoHeader[channel].network);
00786                         // sampling rate
00787                         fsdh.sampleRateFactor=ToNBOs(HWConfig.srExp[channel].base);
00788                         fsdh.sampleRateMultiplier=ToNBOs(-HWConfig.srExp[channel].div);
00789                         seq++;
00790                         /*
00791                         PrintError("seq=%#0x, ndata=%d/%d/%d, N=%d/%d\n",
00792                                         seq,ndata,totalndata,infoHeader[channel].npts,(4*(HEADEROFFSET+16*(nsubframe-1))),BLOCKLEN);
00793                                         */
00794                         frame[HEADEROFFSET+1]=first;
00795                         frame[HEADEROFFSET+2]=previous;
00796                         for (i=HEADEROFFSET; i<BLOCKLEN/4; i++)
00797                         {
00798                                 frame[i]=ToNBOl(frame[i]);
00799                         }
00800                         memcpy(frame,&fsdh,sizeof(struct FixedSeedDataHeader_b1000_b1001));
00801                         // write the block
00802                         fwrite(frame,BLOCKLEN,1,fout);
00803                         // reset the block, frame related data
00804                         memset(frame,0,BLOCKLEN);
00805                         // reset the first value
00806                         first=0xffffffff;
00807                         // reset c,k,nword,subframe,ndata
00808                         c=k=0;
00809                         wndx=2;
00810                         nsubframe=0;
00811                         ndata=0;
00812                 // endif
00813                 }
00814         // while the end of the data is not reached
00815         } while (totalndata!=infoHeader[channel].npts);
00816         fclose(fout);
00817         fclose(mseedChannelFile[channel]);
00818         mseedChannelFile[channel]=NULL;
00819         PrintLog("OUTPUT: %s %d samples (%.2fs)\n",outputFileName,infoHeader[channel].npts,
00820                         (float)infoHeader[channel].npts/(float)infoHeader[channel].sr.npts);
00821         if (noTimeWarning==1)
00822         {
00823                 PrintWarning("The time correction could not be interpolated "
00824                                 "for some blocks of this file. "
00825                                 "Make sure that the time "
00826                                 "correction file exists. "
00827                                 "A preliminary pass with -time may solve the "
00828                                 "problem.\n");
00829         }
00830         unlink(mseedChannelFileName[channel]);
00831 }
00832 
00835 void CloseMSeedFile(int channel)
00836 {
00837         switch (seedDataFormat)
00838         {
00839                 case 3: // int32
00840                         CloseMSeedInt32File(channel);
00841                         break;
00842                 case 10: // steim1
00843                         CloseMSeedSteim1File(channel);
00844                         break;
00845         }
00846 }
00847 
00850 void AddMSeedData(int channel)
00851 {
00852         switch (seedDataFormat)
00853         {
00854                 case 3: // int32
00855                         AddMSeedDataInt32(channel);
00856                         break;
00857                 case 10: // steim1
00858                         AddMSeedDataSteim1(channel);
00859                         break;
00860         }
00861 
00862         if (HWConfig.digitizer.ok>0)
00863         {
00864                 // correct from the ADC delay
00865                 staticDt[channel]=ADCDelay(&(HWConfig.digitizer),channel);
00866 
00867                 // correct from the filter delay
00868                 staticDt[channel]+=FilterDelayTime(&(HWConfig.digitizer),
00869                                 &(HWConfig.srExp[channel]),
00870                                 channel);
00871         }
00872 }
00873