/****************************************************************************/ /* */ /* IAGA_K_index.c */ /* */ /****************************************************************************/ /****************************************************************************/ /* This program computes K indices by the FMI-method for the specified */ /* stations in an IAGA format data file. */ /* */ /* Usage: */ /* IAGA_K_index [-s] [-e] [-o] [-k] [-c] [-a] [-b] [-d] [<] IAGA_file */ /* [-s YYYYMMDD] First day to be computed. If missing then */ /* K indices will be computed all days in the */ /* for which data exists except the first and */ /* the last day. */ /* [-e YYYYMMDD] Last day to be computed. If missing then */ /* only day specified by -s will be computed. */ /* [-o 'Station list'] List of stations delimited by aposthropes */ /* If missing then all stations included. */ /* Stations are identified by three-letter */ /* code and separated by exactly one space. */ /* [-k 'K9 limit list'] List of K=9 limits (in nT) for stations */ /* specified with -o option. Order of stations */ /* must be the same as in -o option. Unit = nT.*/ /* If missing then default values are used. */ /* [-c] Components included in the computation of */ /* the K index. Possible values are */ /* "XY", "X", "Y", "HD", "H" and "D". */ /* If missing then "XY" is used. */ /* Note that it is assumed that the IAGA file */ /* contains components XYZ. H and D are */ /* computed from them. */ /* [-a] Compute the Ak index. If missing then Ak is */ /* not computed. */ /* [-b] Compute the sum of K indices. If missing */ /* then the sum is not computed. */ /* [-d Period] Before K index computation average the data */ /* over period 'Period'. If -d is not defined */ /* then default value is 60 seconds. */ /* IAGA_file Name of IAGA-format file. */ /* */ /****************************************************************************/ /****************************************************************************/ /* Lasse Hakkinen */ /* Finnish Meteorological Institute */ /* Department of Geophysics */ /* P.O.Box 503 */ /* FIN-00101, Helsinki, Finland */ /* e-mail: Lasse.Hakkinen@fmi.fi */ /* phone : (+358)-9-19294634 */ /* fax : (+358)-9-19294603 */ /* */ /* version 1.02 10.01.2020 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.02 10.01.2020 Date strings must be YYYYMMDD (four digits for year) */ /* Replaced StrToSecs -> StrToSecsC and */ /* SecsToStr -> SecsToStrC. */ /* Fixed a bug in handling last day computation. */ /* 1.01 17.10.2017 Added -d option */ /* 1.0 15.06.2010 First official release */ /****************************************************************************/ #include #include #include #include "Usage.h" #include "MagnData.h" #include "K_index.h" #include "IAGA.h" #define OneDay 86400 char *version = "1.02"; char *date = "10.01.2020"; #define HeaderLineCount 5 #define UsageLineCount 22 char *HeaderText[HeaderLineCount] = { "**************************************************************************", " This program computes K indices by the FMI-method for the stations in the", " given IAGA format data file. The file must contain data for the specified", " day and also the previous and following day. ", "**************************************************************************", }; char *UsageText[UsageLineCount] = { " [-s] [-e] [-o] [-k] [-c] [-a] [-b] [<] IAGA_file", " [-s YYYYMMDD] First day to be computed. If missing then K indices", " will be computed for all days for which data exist ", " in the data file, except the first and last day. ", " [-e YYYYMMDD] Last day to be computed. If missing then ", " only day specified with -s option is computed. ", " [-o StationList] List of stations to be included (e.g 'SOR MAS PEL')", " If missing then all stations included. ", " [-k K9 limits] List of K=9 limits (in nT) for the specified ", " stations. Station order must be same as with ", " -o option. ", " [-c] Components Components included in the computation of the K ", " index. Possible values are XY, X, Y, HD, H and D. ", " If missing then XY is used. Note that it is assumed", " that the IAGA file contains components XYZ. If H or", " D is needed then they are computed from XYZ. ", " [-a] Compute the Ak index. ", " [-b] Compute the sum of K indices ", " [-d Period] Before K index computation average the data over ", " period 'Period'. If -d is not defined then default ", " value is 60 seconds. ", " IAGA_file Name of IAGA-format file. " }; /*--------------------------------------------------------------------------*/ /* Find the K=9 limit for the given station. If the StationList is empty */ /* then get the value from the StationInfoTable (defined in StatInfo.h) */ /* otherwise extract the value from the K9List. */ /*--------------------------------------------------------------------------*/ long GetK9Limit(char *StatID,char *StationList, char *K9List) { char *p,*q; /* Dummy pointers */ long found; long value = 0; StationInfoPtr SInfo; if (K9List[0] == '\0') { /* Use the default value */ SInfo = FindStationInfo(StatID); return (10*(SInfo->K9Limit)); } else { /* Use the StationList and K9List */ q = StationList; p = K9List; do { /* Go through all 3-letter ID's */ found = !strncmp(StatID,q,3); if (!found) { while (*p != ' ') p++; while (*p == ' ') p++; } q += 3; } while ((*q++ == ' ') && (!found)); while ((*p != ' ') && (*p != '\0')) { value = 10*value+(long) *p - 48; p++; } return (10*value); } } /*--------------------------------------------------------------------------*/ /* The main procedure */ /*--------------------------------------------------------------------------*/ int main(int argc, char *argv[]) { long params; long FileCount = 0; long status = 0; long j; /* Dummy index */ long K9limit; /* The K=9 limit */ char StartTimeStr[20] = ""; char EndTimeStr[20] = ""; long AveLength = 0; /* Averaging period */ Time_sec StartTime; /* Time of first data block (seconds) */ Time_sec EndTime; /* Time of last block (seconds) */ Time_sec T; /* Current time (seconds) */ char IAGAFileName[200] = ""; char StationList[200] = ""; /* List of stations from the command line */ char K9List[200] = ""; /* List of K=9 limits from the command line */ char TimeStr[20]; /* Time converted to string */ long K_table[8]; /* Table of K indices */ long Ak; /* The Ak number */ long SumK; /* Daily sum of K indices */ long DoAk = 0; /* Flag for computing ak index. */ long DoSumK = 0; /* Flag for computing sum of K indices */ Network_struct NETWORK; StationPtr Station; char Components[10] = "XY"; /* Components included in the computatation */ /*==========================*/ /* 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] != '-') { strcpy(IAGAFileName,argv[params]); FileCount++; } else { switch (*(argv[params]+1)) { case 's' : strcpy(StartTimeStr,argv[++params]); break; case 'e' : strcpy(EndTimeStr,argv[++params]); break; case 'o' : strcpy(StationList,argv[++params]); break; case 'k' : strcpy(K9List,argv[++params]); break; case 'a' : DoAk = 1; break; case 'b' : DoSumK = 1; break; case 'd' : AveLength = atol(argv[++params]); break; case 'c' : strcpy(Components,argv[++params]); 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 the values of the parameters */ /*====================================*/ if (FileCount > 1) { PrintUsage(argv[0],1,HeaderLineCount,UsageLineCount); return FAIL; } /*====================================================*/ /* Read the data from an IAGA format file into memory */ /*====================================================*/ status = ReadIAGA(IAGAFileName,&NETWORK,"","",StationList); if (status != 0) { if (status == FileError) fprintf(stderr,"Error in reading file %s\n",IAGAFileName); if (status == OutOfMemory) fprintf(stderr,"Out of memory while reading IAGA file %s\n",IAGAFileName); if (status == FileFormatError) fprintf(stderr,"%s is not an IAGA format file\n",IAGAFileName); return FAIL; } Station = NETWORK.StationList; if (AveLength == 0) { AveLength = Station->TimeStep; } else { while (Station != NULL) { ComputeStationAverages(Station, AveLength, 10); Station = Station->Next; } } if (strlen(StartTimeStr) == 0) { Station = NETWORK.StationList; SecsToStrC(Station->StartTime + OneDay,StartTimeStr); SecsToStrC(Station->EndTime - 2*OneDay,EndTimeStr); } else if (strlen(EndTimeStr) == 0) strcpy(EndTimeStr,StartTimeStr); /*=======================================*/ /* Go through all the days */ /*=======================================*/ StartTime = StrToSecsC(StartTimeStr,0); EndTime = StrToSecsC(EndTimeStr,0); Station = NETWORK.StationList; if ((StartTime < Station->StartTime + OneDay) || (EndTime > Station->EndTime - 2*OneDay)) { fprintf(stderr,"IAGA file %s does not contain data of previous and next days\n",IAGAFileName); return FAIL; } for (T=StartTime; T<=EndTime; T += OneDay) { /*** Go through all the stations in the IAGA structure ***/ Station = NETWORK.StationList; if ((strcmp(Components,"HD") == 0) || (strcmp(Components,"H") == 0) || (strcmp(Components,"D") == 0)) Compute_HDF(Station); while (Station != NULL) { K9limit = GetK9Limit(Station->StationID,StationList,K9List); if (strcmp(Components,"XY") == 0) K_index(GetDataPtr(Station,T-OneDay,'X'),GetDataPtr(Station,T-OneDay,'Y'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); else if (strcmp(Components,"X") == 0) K_index(GetDataPtr(Station,T-OneDay,'X'),GetDataPtr(Station,T-OneDay,'X'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); else if (strcmp(Components,"Y") == 0) K_index(GetDataPtr(Station,T-OneDay,'Y'),GetDataPtr(Station,T-OneDay,'Y'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); else if (strcmp(Components,"HD") == 0) K_index(GetDataPtr(Station,T-OneDay,'H'),GetDataPtr(Station,T-OneDay,'D'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); else if (strcmp(Components,"H") == 0) K_index(GetDataPtr(Station,T-OneDay,'H'),GetDataPtr(Station,T-OneDay,'H'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); else K_index(GetDataPtr(Station,T-OneDay,'D'),GetDataPtr(Station,T-OneDay,'D'),Station->TimeStep,K9limit,(Station->Longitude)/100,MissingValue,K_table); /* print the results */ SecsToStrC(T,TimeStr); TimeStr[8] = '\0'; printf("%s %s ",TimeStr,Station->StationID); for (j=0;j<8;j++) { if (K_table[j] < 0) printf(" *"); else printf("%2d",K_table[j]); } if (DoSumK == 1) { SumK = ComputeSumK(K_table); if (SumK < 0) printf(" *"); else printf("%5d",SumK); } if (DoAk == 1) { Ak = ComputeAk(K_table); if (Ak < 0) printf(" *"); else printf("%5d",Ak); } printf("\n"); Station = Station->Next; } } FreeNetwork(&NETWORK); return OK; }