Titan2Reader: trigger.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 <stdlib.h>
00024 #include <stdio.h>
00025 #include <time.h>
00026 #include <ctype.h>
00027 #include <math.h>
00028 #include <errno.h>
00029 #include <unistd.h>
00030 #include <string.h>
00031 #include <sys/types.h>
00032 #include <sys/stat.h>
00033 #include <fcntl.h>
00034 #include <dirent.h>
00035 #include <assert.h>
00036 #include <signal.h>
00037 
00038 #include "common.h"
00039 
00040 int firstevt=0;
00041 
00042 float sta=1., lta=60., hp=0.05, lp=20.0, threshold=2.0, triggerMinLen=-1.;
00043 
00044 char outfile[128] = "";
00045 FILE *fout;
00046 
00047 char evtfile[128] = "";
00048 FILE *fevt;
00049 
00050 
00051 struct trigger_t {
00052         struct Titan2TimeType start;
00053         float ksta,klta,kHP,kLP;
00054         float LP;
00055         float HP;
00056         float sta0;
00057         float lta0;
00058         int64_t offset0;
00059         float sta;
00060         float lta;
00061 } trigger[MAX_CHANNEL];
00062 
00063 int ParseTriggerOptions(char *cmd)
00064 {
00065         char *s=cmd;
00066         int next,error,i;
00067         for (i=0; i<MAX_CHANNEL; i++)
00068         {
00069                 memset(&trigger[i],0,sizeof(struct trigger_t));
00070         }
00071         if (*s==0)
00072         {
00073                 return(0);
00074         }
00075         if (*s!=',')
00076         {
00077                 return(-1);
00078         }
00079         while (*s)
00080         {
00081                 PrintDebug("%s\n",s);
00082                 s++;
00083                 next=0;
00084                 error=__LINE__;
00085                 if (strncmp(s,"lp=",3)==0)
00086                 {
00087                         error=__LINE__;
00088                         if (sscanf(s+3,"%f%n",&lp,&next)!=1)
00089                         {
00090                                 break;
00091                         }
00092                         error=0;
00093                         next+=3;
00094                 }
00095                 if (strncmp(s,"hp=",3)==0)
00096                 {
00097                         error=__LINE__;
00098                         if (sscanf(s+3,"%f%n",&hp,&next)!=1)
00099                         {
00100                                 break;
00101                         }
00102                         error=0;
00103                         next+=3;
00104                         if (hp<0.)
00105                         {
00106                                 hp=0.002;
00107                         }
00108                 }
00109                 if (strncmp(s,"m=",2)==0)
00110                 {
00111                         error=__LINE__;
00112                         if (sscanf(s+2,"%f%n",&triggerMinLen,&next)!=1)
00113                         {
00114                                 break;
00115                         }
00116                         error=0;
00117                         next+=2;
00118                 }
00119                 if (strncmp(s,"t=",2)==0)
00120                 {
00121                         error=__LINE__;
00122                         if (sscanf(s+2,"%f%n",&threshold,&next)!=1)
00123                         {
00124                                 break;
00125                         }
00126                         error=0;
00127                         next+=2;
00128                 }
00129                 if (strncmp(s,"l=",2)==0)
00130                 {
00131                         error=__LINE__;
00132                         if (sscanf(s+2,"%f%n",&lta,&next)!=1)
00133                         {
00134                                 break;
00135                         }
00136                         error=0;
00137                         next+=2;
00138                 }
00139                 if (strncmp(s,"s=",2)==0)
00140                 {
00141                         error=__LINE__;
00142                         if (sscanf(s+2,"%f%n",&sta,&next)!=1)
00143                         {
00144                                 break;
00145                         }
00146                         error=0;
00147                         next+=2;
00148                 }
00149                 if (strncmp(s,"e=",2)==0)
00150                 {
00151                         int i=0;
00152                         error=__LINE__;
00153                         s+=2;
00154                         while ((*s!=0)&&(*s!=','))
00155                         {
00156                                 evtfile[i++]=*s;
00157                                 s++;
00158                         }
00159                         evtfile[i]=0;
00160                         error=0;
00161                 }
00162                 if (strncmp(s,"o=",2)==0)
00163                 {
00164                         int i=0;
00165                         error=__LINE__;
00166                         s+=2;
00167                         while ((*s!=0)&&(*s!=','))
00168                         {
00169                                 outfile[i++]=*s;
00170                                 s++;
00171                         }
00172                         outfile[i]=0;
00173                         error=0;
00174                 }
00175                 if (error)
00176                 {
00177                         PrintError("could not parse trigger option string (error %d)\n",error);
00178                         return(-1);
00179                 }
00180                 s+=next;
00181         }
00182         fout=stdout;
00183         fevt=NULL;
00184         if (strlen(evtfile)>0)
00185         {
00186                 fevt=fopen(evtfile,"at");
00187                 if (fevt==NULL)
00188                 {
00189                         PrintError("%s: %m\n",evtfile);
00190                 }
00191         }
00192         if (strlen(outfile)>0)
00193         {
00194                 fout=fopen(outfile,"at");
00195                 if (fout==NULL)
00196                 {
00197                         PrintError("%s: %m\n",outfile);
00198                         fout=stdout;
00199                 }
00200         }
00201         firstevt=1;
00202         return(0);
00203 }
00204 
00205 
00206 
00207 void TriggerOneSample(channel, value)
00208 {
00209         float F;
00210         // pour chaque échantillon Y des canaux actifs en trigger:
00211         // LP=LP+(Y-LP)*KLP
00212         if (lp>0.)
00213         {
00214                 trigger[channel].LP+=((float)value-trigger[channel].LP)*
00215                         trigger[channel].kLP;
00216         }
00217         else
00218         {
00219                 trigger[channel].LP=(float)value;
00220         }
00221         // HP=HP+(Y-HP)*KHP
00222         if (hp>0.)
00223         {
00224                 trigger[channel].HP+=((float)value-trigger[channel].HP)*
00225                         trigger[channel].kHP;
00226         }
00227         else
00228         {
00229                 trigger[channel].HP=0.0;
00230         }
00231         // F=LP-HP
00232         F=trigger[channel].LP-trigger[channel].HP;
00233         //printf("%f %d\n",F,value);
00234         // sta=sta+(abs(F)-sta)*ksta
00235         trigger[channel].sta+=(fabs(F)-trigger[channel].sta)*
00236                 trigger[channel].ksta;
00237         // lta=lta+(abs(F)-lta)*klta
00238         trigger[channel].lta+=(fabs(F)-trigger[channel].lta)*
00239                 trigger[channel].klta;
00240         // si canal déjà déclenché
00241         if ((trigger[channel].start.usecond>0)&&
00242                         ((correctedTime.usecond-trigger[channel].start.usecond)>>32>0))
00243         {
00244                 // traiter le sta0 ???
00245                 // si sta/lta>seuil
00246                 if (trigger[channel].sta>threshold*trigger[channel].lta)
00247                 {
00248                         // rien
00249                 }
00250                 else
00251                 // sinon
00252                 {
00253                         int duration=(correctedTime.usecond-trigger[channel].start.usecond)>>32;
00254                         // si assez long
00255                         if (duration>=triggerMinLen)
00256                         {
00257                                 // pousser le trigger pour ce canal
00258                                 if (fevt!=NULL)
00259                                 {
00260                                         struct Titan2TimeType t;
00261                                         memcpy(&t,&trigger[channel].start,sizeof(struct Titan2TimeType));
00262                                         t.usecond-=(30LL<<32);
00263                                         fprintf(fevt,"%s.000 %d\n",StrTimeNS(&t),180);
00264                                         fflush(fevt);
00265                                 }
00266                                 if (firstevt==1)
00267                                 {
00268                                         fprintf(fout,"#TRIGGER station "
00269                                                      "S/N ch date/time "
00270                                                      "epoch duration sta "
00271                                                      "lta fileoffset\n");
00272                                 }
00273                                 firstevt=0;
00274                                 fprintf(fout,"TRIGGER %-9s %-6s %2d %s %8lld %3d %.0f %.0f %lld\n",
00275                                                 HWConfig.fieldIdent,
00276                                                 HWConfig.serialNumber,
00277                                                 channel,
00278                                                 StrTimeNS(&trigger[channel].start),
00279                                                 trigger[channel].start.usecond>>32,
00280                                                 duration,
00281                                                 trigger[channel].sta0,
00282                                                 trigger[channel].lta0,
00283                                                 trigger[channel].offset0);
00284                                 // reinit du trigger pour ce canal
00285                         }
00286                         // finsi
00287                         trigger[channel].start.usecond=0;
00288                 // finsi
00289                 }
00290         }
00291         // sinon
00292         else
00293         {
00294                 // si sta/lta>seuil
00295                 if (trigger[channel].sta>threshold*trigger[channel].lta)
00296                 {
00297                         // déclencher le canal
00298                         trigger[channel].start.usecond=correctedTime.usecond;
00299                         trigger[channel].sta0=trigger[channel].sta;
00300                         trigger[channel].lta0=trigger[channel].lta;
00301                         trigger[channel].offset0=inputOffset;
00302                 // finsi
00303                 }
00304         // finsi
00305         }
00306 }
00307 
00308 float FilterT2K(int channel,float T)
00309 {
00310         float dummy=2.;
00311         dummy=T*(float)HWConfig.srExp[channel].base/
00312                 (float)HWConfig.srExp[channel].div;
00313         return(1.-exp(-1./dummy));
00314 }
00315 
00316 // format de sortie:
00317 // STATION SN CANAL YYYY.MM.DD.hh.mm.ss.sss epoch.sss durée sta lta
00318 
00319 void AddTriggerData(int channel)
00320 {
00321         int i;
00322         if ((HWConfig.srObs[channel].npts==0)||
00323                         (HWConfig.srExp[channel].npts==0))
00324         {
00325                 return;
00326         }
00327         trigger[channel].ksta=FilterT2K(channel,sta);
00328         trigger[channel].klta=FilterT2K(channel,lta);
00329         trigger[channel].kHP=FilterT2K(channel,1./hp);
00330         trigger[channel].kLP=FilterT2K(channel,1./lp);
00331 
00332         for (i=0; i<HWConfig.srObs[channel].npts; i++)
00333         {
00334                 TriggerOneSample(channel, HWConfig.data[channel][i]);
00335         }
00336 }
00337