/****************************************************************************/ /* */ /* IAGA_time_check.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program reads two IAGA format data files and compares the timing */ /* between the files by cross-correlation techniques. */ /* */ /* Usage: */ /* IAGA_time_check -h -c File1 File2 */ /* -c Components Components to be compared. */ /* -h HourCount Number of hours in used in a single comparison */ /* Allowed values are 1,2,3,4,6,12 or 24. */ /* File1 Name of the first IAGA-format file. Gaps in */ /* this file are filled with the data values from */ /* the second file. */ /* File2 Name of the second IAGA-format file. */ /* NewIAGAFile Name of combined IAGA-format file. */ /* */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Geophysical Research Division */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.03 13.02.2020 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.03 13.02.2020 Date strings must be YYYYMMDD (four digits for year) */ /* Replaced StrToSecs -> StrToSecsC and */ /* SecsToStr -> SecsToStrC. */ /* 1.02 09.09.1999 Fixed a Y2K bug in NewTime.h file which resulted in */ /* incorrect year in date strings if year >= 2000. */ /* 1.01 12.11.1997 Code polished. */ /* 1.0 05.11.1997 First official release */ /****************************************************************************/ #include #include #include #include #include "Usage.h" #include "MagnData.h" #include "IAGA.h" char *version = "1.03"; char *date = "13.02.2020"; #define HeaderLineCount 4 #define UsageLineCount 10 #define MaxTimeOffset 1200 char *HeaderText[HeaderLineCount] = { "**************************************************************************", " IAGA_time_check reads in two IAGA format data files and compares the ", " timing between the two files. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " -c -h File1 File2 ", " -c ComponentList Timing comparison is made separately for each ", " specified components. Default value is 'XYZ'. ", " -h HourCount Data is divided into HourCount long units. ", " Timing comparison is made separately for each of ", " these subintervals. Possible values for HourCount", " are 1,2,3,4,6,12 or 24. Default is 2. ", " File1 Name of first IAGA-format file. Data gaps in this", " file are filled with the values from file #2. ", " File2 Name of the second IAGA-format file. ", }; /*--------------------------------------------------------------------------*/ /* Compute correlation for one component and for 'HourCount' hours */ /* starting from hour 'Hour'. */ /*--------------------------------------------------------------------------*/ long ComputeCorrelation(StationPtr S1, StationPtr S2, char Comp, long Hour, long HourCount) { long *DataPtr1,*DataPtr2; double *CurrPtr; double *DataArea; long Count, Corr_count; long F1,F2; double Sum, Diff; Time_sec T1,T2,BaseTime, tMin; long dT; long DataSize,i,j; long valid; double Min = DBL_MAX; long B1,B2; /* Allocate space for the correlation vector */ Corr_count = 2*MaxTimeOffset/(S1->TimeStep) + 1; DataSize = (sizeof(double))*Corr_count; DataArea = (double *) malloc(DataSize); CurrPtr = DataArea; // Count = 3600*HourCount/S1->TimeStep; Count = (3600*HourCount - 2*MaxTimeOffset)/S1->TimeStep; BaseTime = S1->StartTime + 3600*Hour + MaxTimeOffset; B1 = ComputeAverage_Conditional(S1, Comp, BaseTime, 3600*HourCount, 50); B2 = ComputeAverage_Conditional(S2, Comp, BaseTime, 3600*HourCount, 50); if ((B1 == MissingValue) || (B2 == MissingValue)) return (MissingValue); /* Compute the correlation vector */ for (dT = -MaxTimeOffset; dT <= MaxTimeOffset; dT += S1->TimeStep) { Sum = 0.0; valid = 0; T1 = BaseTime; T2 = BaseTime + dT; for (i = 0; i < Count; i++) { F1 = GetDataValue(S1, T1, Comp); F2 = GetDataValue(S2, T2, Comp); if ((F1 != MissingValue) && (F2 != MissingValue)) { Diff = (double)((F1 - B1) - (F2 - B2)); Sum += Diff*Diff; valid++; } else { } T1 += S1->TimeStep; T2 += S1->TimeStep; } if (valid < Count/2) return (MissingValue); Sum = Sum/valid; *CurrPtr = Sum; CurrPtr++; if (Sum < Min) { Min = Sum; tMin = dT; } } free(DataArea); return (tMin); } /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; long status = 0; long FileCount = 0; char FileName[100] = ""; /* Name of the original data file */ char CompStr[5] = "XYZ"; /* List of components */ char Component; Network_struct IMAGE1; /* Structure containing the data */ /* for station in the first file. */ Network_struct IMAGE2; /* Structure containing the data */ /* for station in the second file. */ StationPtr Station1; /* Pointer to the first station */ StationPtr Station2; /* Pointer to the second station */ long Comp; /* Current magnetic component */ long Shift; char DateStr[20]; char YearStr[3]; char MonthStr[3]; char DayStr[3]; long Hour, EndHour; long HourCount = 2; /* Number of hours used in one comparison */ long PrintFormat = 0; /*==========================*/ /* Analyse the command line */ /*==========================*/ if (argc == 1) { /* No arguments, show the usage text */ PrintUsage(argv[0],0,HeaderLineCount,UsageLineCount); return OK; } for (params = 1; params < argc; params++) { if (*argv[params] != '-') { argv[FileCount++] = argv[params]; } else { switch (*(argv[params]+1)) { case 'c' : strcpy(CompStr,argv[++params]); break; case 'h' : HourCount = atol(argv[++params]); break; case 'a' : PrintFormat = 1; break; default : fprintf(stderr,"\n### %s: \"%s\" is not an option.\n\n", argv[0], argv[params]); PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; break; } } } /*=================================================*/ /* Check that all necessary parameters are defined */ /*=================================================*/ if (FileCount != 2) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } /*======================================================*/ /* Read data from the two IAGA format files into memory */ /*======================================================*/ strcpy(FileName,argv[0]); status = ReadIAGA(FileName,&IMAGE1,NULL,NULL,NULL); if (status != 0) { if (status == FileError) fprintf(stderr,"Error in reading file %s\n",FileName); if (status == OutOfMemory) fprintf(stderr,"Out of memory while reading file %s\n",FileName); if (status == FileFormatError) fprintf(stderr,"%s is not an IAGA format file\n",FileName); return FAIL; } strcpy(FileName,argv[1]); status = ReadIAGA(FileName,&IMAGE2,NULL,NULL,NULL); if (status != 0) { if (status == FileError) fprintf(stderr,"Error in reading file %s\n",FileName); if (status == OutOfMemory) fprintf(stderr,"Out of memory while reading file %s\n",FileName); if (status == FileFormatError) fprintf(stderr,"%s is not an IAGA format file\n",FileName); return FAIL; } /*=======================================================*/ /* Check that the values of the parameters are resonable */ /*=======================================================*/ if ((HourCount != 1) && (HourCount != 2) && (HourCount != 3) && (HourCount != 4) && (HourCount != 6) && (HourCount != 12) && (HourCount != 24)) { fprintf(stderr,"### Illegal hour count option -h : %d\n",HourCount); return FAIL; } /*======================================*/ /* Write header */ /*======================================*/ Station1 = IMAGE1.StationList; Station2 = IMAGE2.StationList; SecsToStrC(Station1->StartTime,DateStr); DateStr[8] = '\0'; if (PrintFormat == 0) { printf("%s - %s %s\n",Station1->StationID,Station2->StationID,DateStr); } /*======================================*/ /* Go through all specified components */ /*======================================*/ EndHour = (Station2->EndTime - Station2->StartTime)/3600; SecsToStrC(Station1->StartTime, DateStr); strncpy(YearStr, DateStr,4); strncpy(MonthStr, DateStr+4,2); strncpy(DayStr, DateStr+6,2); YearStr[4] = '\0'; MonthStr[2] = '\0'; DayStr[2] = '\0'; for (Comp=0;Comp MaxTimeOffset) if (PrintFormat == 1) printf(" NaN"); else printf(" ---",Shift); else if (PrintFormat == 1) printf("%5d", Shift); else printf(" %4d",Shift); } if (PrintFormat == 1) printf("\n"); } if (PrintFormat == 0) printf("\n"); } /* End of Component Loop */ FreeNetwork(&IMAGE1); FreeNetwork(&IMAGE2); return OK; }