Code:
#define MAIN
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "/media/Elements/ANT/progs_HLP/HEADFILE/mysac.h"
#include "string.h"
/* a program to take the list files, for up to 30 files,
and stack all the files that are present in all
directories. The list files must have the absolute
path as the first line. The other lines should be
a file name followed by a tab then the number of days
of raw data that lead to this correlation. */
/* information about all files in all directories are
read into this af[] structure. */
// change the delta of sac file when it is not 1.0
/* _cv vision dosn't use the SAC program but use the pure C progrrame.
*
*/
struct onefile { /* file name and number */
char name[100]; /* complete file name */
int ndays; /* number of days for the given correlation */
int onoff; /* tells if file has been stacked (1) or not (0) */
float depmin; /* min signal value, from sac header depmin value */
float depmax; /* like depmin, max signal value */
};
struct allfiles { /*for whole filelist with dir and len*/
char dir[100]; /* absolute path, ending in a "/" */
struct onefile fname[400000]; /* structure for each file */
int len; /* total number of files in a given directory */
} af[40]; /* up to 30 month.lst files can be used, this can be changed */
/**/
void write_sac (char *fname, float *sig, SAC_HD *SHD)
/*----------------------------------------------------------------------------
* ----------------------------------------------------------------------------*/
{
FILE *fsac;
int i;
/*..........................................................................*/
fsac = fopen(fname, "wb");
if ( !SHD ) SHD = &SAC_HEADER;
SHD->iftype = (long)ITIME;
SHD->leven = (long)TRUE;
SHD->lovrok = (long)TRUE;
SHD->internal4 = 6L;
/*+++++++++++++++++++++++++++++++++++++++++*/
SHD->depmin = sig[0];
SHD->depmax = sig[0];
for ( i = 0; i < SHD->npts ; i++ )
{
if ( SHD->depmin > sig[i] ) SHD->depmin = sig[i];
if ( SHD->depmax < sig[i] ) SHD->depmax = sig[i];
}
fwrite(SHD,sizeof(SAC_HD),1,fsac);
fwrite(sig,sizeof(float),(int)(SHD->npts),fsac);
fclose (fsac);
}
/**/
/*c/////////////////////////////////////////////////////////////////////////*/
/*--------------------------------------------------------------------------*/
SAC_HD *read_sac (char *fname, float *sig, SAC_HD *SHD, int nmax)
/*----------------------------------------------------------------------------
----------------------------------------------------------------------------*/
{
FILE *fsac;
/*..........................................................................*/
if((fsac = fopen(fname, "rb")) == NULL) {
printf("could not open sac file to read%s \n", fname);
exit(1);
}
if ( !fsac )
{
/*fprintf(stderr,"file %s not find\n", fname);*/
return NULL;
}
if ( !SHD ) SHD = &SAC_HEADER;
fread(SHD,sizeof(SAC_HD),1,fsac);
if ( SHD->npts > nmax )
{
fprintf(stderr,
"ATTENTION !!! dans le fichier %s npts est limite a %d",fname,nmax);
SHD->npts = nmax;
}
fread(sig,sizeof(float),(int)(SHD->npts),fsac);
fclose (fsac);
/*------------- calcule de t0 ----------------*/
{
int eh, em ,i;
float fes;
char koo[9];
for ( i = 0; i < 8; i++ ) koo[i] = SHD->ko[i];
koo[8] = 0 ;
SHD->o = SHD->b + SHD->nzhour*3600. + SHD->nzmin*60 +
SHD->nzsec + SHD->nzmsec*.001;
sscanf(koo,"%d%*[^0123456789]%d%*[^.0123456789]%g",&eh,&em,&fes);
SHD->o -= (eh*3600. + em*60. + fes);
/*-------------------------------------------*/}
return SHD;
}
#define SLEN 400000
float sig0[SLEN];
float cvsig[SLEN];
float tempsig[SLEN];
int main(int argc, char *argv[])
{
FILE *fp1, *fp0;
char ch;
int i = 0, j = 0, k = 0, quit = 0, nfiles = 0;
int inotbase = 0, jbase = 0, jnotbs = 0;
int div_by = 0, mindays = 0, ibase = 0;
char dir[100], filename[100], nombre[100];
char dump[100], outdir[100], fname[100];
float minv = 2000000;
SAC_HD cvSAC_HEADER, tempSAC_HEADER;
/* check for comand line arg */
if(argc!=4) {
printf("Usage: enter filelist outfile_directory mindays\n");
exit(1);
}
strcpy(outdir, argv[2]);
/* open file for reading */
if((fp0 = fopen(argv[1], "r"))==NULL) {
printf("cannot open infile.\n");
exit(1);
}
fgets(filename, 60, fp0); /*name of month.lst */
/* fill the structure with all the file information for up to
30 files */
do{
i = strlen(filename);
filename[i-1] ='\0';
printf("file name %d is %s\n", k, filename);
/* open file for reading */
if((fp1 = fopen(filename, "r"))==NULL) {
printf("cannot open infile number %d.\n", k);
exit(1);
}
/* first line of file must be the absolute path*/
fgets(dir, 80, fp1);
i = strlen(dir);
dir[i-1] = '\0';
strcpy(af[k].dir, dir);
//fprintf(stderr,"the SAC files are at : %s \n",dir);
if( fgets(nombre,80,fp1)==NULL)
goto next;
do {
i = strlen(nombre);
nombre[i-1] = '\0';
sprintf(fname, "%s%s", af[k].dir, nombre);
if (fmod(j,50.0)==0.0)
//fprintf(stderr,"%s%s\n",af[k].dir,nombre);
/*---------------- reading sac file -------------------*/
if ( read_sac (fname, sig0, &SAC_HEADER, SLEN) == NULL ){
fprintf(stderr,"file %s not found\n", fname);
return 0;
}
/* judge if the sac file is a null file with only headfile no signal */
if( SAC_HEADER.depmin < minv ) {
strcpy(af[k].fname[j].name, nombre);
/* get the number of days */
// fgets(dump, 30, fp1);
//af[k].fname[j].ndays = atoi(dump);
af[k].fname[j].ndays = 20;
af[k].fname[j].onoff = 0;
af[k].fname[j].depmin = ((-1)*SAC_HEADER.depmin);
af[k].fname[j].depmax = SAC_HEADER.depmax;
j++;
}
else {
//fgets(dump, 30, fp1);
}
fgets(nombre,80,fp1);
} while(!feof(fp1));
next:
/* j is number of files for a given dir */
af[k].len = j;
fprintf(stderr, "length of file %d is %d\n", k, af[k].len);
j = 0;
k++;
fgets(filename, 60, fp0);
} while(!feof(fp0));
fprintf(stderr, "all info read fine\n\n");
nfiles = k; /* number of files as determined from loop*/
fclose(fp1);
fclose(fp0);
/*next section checks for file existence in all directories
and writes script to add the files together. It will only be
executed if the file is in all directories. */
/* inotbase is file num. inotbase=0 is reference file*/
jbase = 0;/* j is line of first (ibase) file */
jnotbs = 0;/* jnotbs is line of second file */
char cvsacname[30];
char tempsacname[30];
char cvobsac[30];
int cvi=0;
/* increment over nfiles so if a file is in the second, third
etc but not the first it will not be missed */
for(ibase = 0; ibase < nfiles; ibase++){
if (af[ibase].len==0)
continue;
do{ /* outer do while */
inotbase = ibase + 1;
if(af[ibase].fname[jbase].onoff==0){
/* write sac script header */
/*if((fp1 = fopen("dostack.csh", "w"))==NULL) {
fprintf(stderr, "cannot open dostack.csh \n");
exit(1);
}
fprintf(stderr, "filename is %s\n", af[ibase].fname);*/
// read sac file to cvsig[200000], SAC_HD cvsac_hd;
sprintf(cvsacname,"%s%s\0",af[ibase].dir, af[ibase].fname[jbase].name);
// fprintf(stderr,"%s\n",cvsacname);
if ( read_sac (cvsacname, cvsig, &cvSAC_HEADER, SLEN) == NULL ){
fprintf(stderr,"file %s not found\n", af[ibase].fname);
return 0;
}
if (cvSAC_HEADER.delta!=1.0) cvSAC_HEADER.delta=1.0;
af[ibase].fname[jbase].onoff=1;
//fprintf(fp1, "#!/bin/csh\n");
//fprintf(fp1, "sac << END\n");
//fprintf(fp1, "r %s%s\n", af[ibase].dir, af[ibase].fname[jbase].name);
// if(af[ibase].fname[jbase].depmin > af[ibase].fname[jbase].depmax){
// fprintf(fp1, "div %f\n", af[ibase].fname[jbase].depmin);
// }
// else fprintf(fp1, "div %f\n", af[ibase].fname[jbase].depmax);
//fprintf(fp1, "mul %d\n", af[ibase].fname[jbase].ndays);
//fprintf(fp1, "w a\n\n");
/* div_by will sum up total number of days from each month */
div_by = af[ibase].fname[jbase].ndays;
do{ /* inner do while */
/* if strings are different, increment jnotbs, quit if jnotbs
gets to reach the end of the file */
if(strcmp(af[ibase].fname[jbase].name, af[inotbase].fname[jnotbs].name)) {
if(jnotbs == af[inotbase].len) {
fprintf(stderr, "match wrong\n");
if(inotbase != nfiles) inotbase++;
jnotbs = 0;
}
else jnotbs++;
}
/* if file names are the same, write this part of the script*/
else {
//fprintf(fp1, "r %s%s\n", af[inotbase].dir, af[ibase].fname[jbase].name);
sprintf(tempsacname,"%s%s\0", af[inotbase].dir, af[ibase].fname[jbase].name);
fprintf(stderr,"%s\n",tempsacname);
if ( read_sac (tempsacname, tempsig, &tempSAC_HEADER, SLEN) == NULL ){
fprintf(stderr,"file %s not found\n", af[ibase].fname);
return 0;
}
if (tempSAC_HEADER.delta!=1.0) tempSAC_HEADER.delta=1.0;
af[inotbase].fname[jnotbs].onoff = 1;
// if(af[inotbase].fname[jnotbs].depmin > af[inotbase].fname[jnotbs].depmax){
// fprintf(fp1, "div %f\n", af[inotbase].fname[jnotbs].depmin);
// }
//else fprintf(fp1, "div %f\n", af[inotbase].fname[jnotbs].depmax);
//fprintf(fp1, "mul %d\n", af[inotbase].fname[jnotbs].ndays);
//fprintf(fp1, "w aa\n");
//fprintf(fp1, "r a\n");
//fprintf(fp1, "addf aa\n");
//fprintf(fp1, "w a\n\n");
fprintf(stderr,"npts: %d\n",tempSAC_HEADER.npts);
if (tempSAC_HEADER.npts == cvSAC_HEADER.npts) {
for (cvi=0; cvi<tempSAC_HEADER.npts; cvi++) {
cvsig[cvi]=cvsig[cvi]+tempsig[cvi];
}
div_by += af[inotbase].fname[jnotbs].ndays;
}
else fprintf(stderr,"the sac file %s is not consistant with sac file %s\n", cvsacname,tempsacname);
inotbase++;
jnotbs = 0;
}
} while(inotbase != nfiles); /* inner do while */
//fprintf(fp1, "r a\n");
//fprintf(fp1, "ch user0 %d\n", div_by);
//fprintf(fp1, "w %s%s\n", outdir, af[ibase].fname[jbase].name);
//fprintf(fp1, "END\n\n");
//fclose(fp1);
/* if inotbase = nfiles, file present in all directories */
/* this could be made into something less than nfiles
for example, if nfiles = 12, this could be executed
for inotbase > 10 or something. */
if(inotbase == nfiles){
//fprintf(stderr, "file %s present in all dirs \n", af[ibase].fname[jbase].name);
}
strcpy(dump, argv[3]);
mindays=atoi(dump);
if(div_by >= mindays){
//system("csh dostack.csh");
//sprintf(dump, "mv dostack.csh scripts/%d%s.csh", ibase, af[ibase].fname[jbase].name);
//system(dump);
//fprintf(stderr,"outdir : %s\n",outdir);
sprintf(cvobsac,"%s%s\0",outdir, af[ibase].fname[jbase].name);
cvSAC_HEADER.delta=1.0;
write_sac (cvobsac,&(cvsig[0]),&cvSAC_HEADER);
}
} /* end if loop */
/* reset the indices */
quit = 0;
jnotbs = 0;
jbase++;
} while(jbase != af[ibase].len); /* outer do while */
jbase = 0;
} /* end for loop */
return 0;
}
The mysac.h file is
Code:
#ifndef MYSAC
#define MYSAC
typedef struct sac {
float delta, depmin, depmax, scale, odelta;
float b, e, o, a, internal1;
float t0, t1, t2, t3, t4;
float t5, t6, t7, t8, t9;
float f, resp0, resp1, resp2, resp3;
float resp4, resp5, resp6, resp7, resp8;
float resp9, stla, stlo, stel, stdp;
float evla, evlo, evel, evdp, unused1;
float user0, user1, user2, user3, user4;
float user5, user6, user7, user8, user9;
float dist, az, baz, gcarc, internal2;
float internal3, depmen, cmpaz, cmpinc, unused2;
float unused3, unused4, unused5, unused6, unused7;
float unused8, unused9, unused10, unused11, unused12;
int nzyear, nzjday, nzhour, nzmin, nzsec;
int nzmsec, internal4, internal5, internal6, npts;
int internal7, internal8, unused13, unused14, unused15;
int iftype, idep, iztype, unused16, iinst;
int istreg, ievreg, ievtyp, iqual, isynth;
int unused17, unused18, unused19, unused20, unused21;
int unused22, unused23, unused24, unused25, unused26;
int leven, lpspol, lovrok, lcalda, unused27;
char kstnm[8], kevnm[16];
char khole[8], ko[8], ka[8];
char kt0[8], kt1[8], kt2[8];
char kt3[8], kt4[8], kt5[8];
char kt6[8], kt7[8], kt8[8];
char kt9[8], kf[8], kuser0[8];
char kuser1[8], kuser2[8], kcmpnm[8];
char knetwk[8], kdatrd[8], kinst[8];
} SAC_HD;
#ifndef MAIN
extern SAC_HD SAC_HEADER;
#else
SAC_HD SAC_HEADER;
#endif
static SAC_HD sac_null = {
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345., -12345., -12345., -12345., -12345.,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
-12345, -12345, -12345, -12345, -12345,
{ '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ',' ',' ',' ',' ',' ',' ',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }, { '-','1','2','3','4','5',' ',' ' },
{ '-','1','2','3','4','5',' ',' ' }
};
/* defines for logical data types */
#define TRUE 1
#define FALSE 0
/* defines for enumerated data types */
#define IREAL 0
#define ITIME 1
#define IRLIM 2
#define IAMPH 3
#define IXY 4
#define IUNKN 5
#define IDISP 6
#define IVEL 7
#define IACC 8
#define IB 9
#define IDAY 10
#define IO 11
#define IA 12
#define IT0 13
#define IT1 14
#define IT2 15
#define IT3 16
#define IT4 17
#define IT5 18
#define IT6 19
#define IT7 20
#define IT8 21
#define IT9 22
#define IRADNV 23
#define ITANNV 24
#define IRADEV 25
#define ITANEV 26
#define INORTH 27
#define IEAST 28
#define IHORZA 29
#define IDOWN 30
#define IUP 31
#define ILLLBB 32
#define IWWSN1 33
#define IWWSN2 34
#define IHGLP 35
#define ISRO 36
#define INUCL 37
#define IPREN 38
#define IPOSTN 39
#define IQUAKE 40
#define IPREQ 41
#define IPOSTQ 42
#define ICHEM 43
#define IOTHER 44
#define IGOOD 45
#define IGLCH 46
#define IDROP 47
#define ILOWSN 48
#define IRLDTA 49
#define IVOLTS 50
#define INIV51 51
#define INIV52 52
#define INIV53 53
#define INIV54 54
#define INIV55 55
#define INIV56 56
#define INIV57 57
#define INIV58 58
#define INIV59 59
#define INIV60 60
#define FCS "%15.7f%15.7f%15.7f%15.7f%15.7f\n"
#define ICS "%10d%10d%10d%10d%10d\n"
#define CCS1 "%-8.8s%-8.8s%-8.8s\n"
#define CCS2 "%-8.8s%-16.16s\n"
SAC_HD *read_sac (char *fname, float *sig, SAC_HD *SHD, int nmax);
void write_sac (char *fname, float *sig, SAC_HD *SHD);
#endif
Originally Posted by Adak
