/****************************************************************************/ /* */ /* MagnData.h */ /* */ /****************************************************************************/ /****************************************************************************/ /* This is a C language header file that defines data structures and */ /* routines for manipulating magnetic data from an arbitrary number of */ /* recording stations. These routines are used in manipulating the data in */ /* memory, separate routines (e.g. for IAGA and GADF format data) must be */ /* used to read the data from files into these data structures. */ /* The principal idea in designing these structures is that they must be */ /* able to handle data from arbitrary number of stations and for arbitrary */ /* time periods. Therefore space for these structures must be allocated */ /* dynamically during program execution and the stations must be defined */ /* as a linked list of data structures. */ /****************************************************************************/ /****************************************************************************/ /* 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.14 09.01.2008 */ /****************************************************************************/ /****************************************************************************/ /* Version history: */ /* */ /* 1.14 09.01.2008 Added ScaleData-routine */ /* 1.13 10.08.2005 Added ComputeAverage_Conditional-routine */ /* 1.12 07.09.2004 Added components H, D and F into Station_struct and */ /* added Compute_HDF -routine. */ /* 1.11 03.08.2004 Added component B into Station_struct. */ /* 1.10 08.11.2000 Added FillWithMissingData-routine and modified */ /* NewOneDayBlock-routine so that it initializes data with */ /* missing data values. */ /* 1.09 14.01.2000 Added RemoveStation-routine. */ /* 1.08 15.01.1999 Added SetDataValue-routine and modified InterpolateData */ /* function. */ /* 1.07 27.10.1998 Added a new component into Station_struct. This is used */ /* e.g. in Dump-files. */ /* 1.06 23.10.1998 Added ChangeUnitUP and ChangeUnitDOWN-routines. */ /* Also changed FindMaxMin so that it returns MissingValue */ /* for Max and Min if all data points are missing */ /* 1.05 27.10.1997 Added NewStation and CopyMagnData-routines */ /* 1.04 20.10.1997 Added InterpolateData-routine */ /* 1.03 7.04.1997 Fixed a bug in FindMaxMin routine which gave incorrect */ /* results if the magnetic field was monotonous (either */ /* growing or diminishing) in the specified interval. */ /* 1.02 19.02.1997 Added FindMaxAmplitude routine */ /* 1.01 9.11.1996 modified ComputeAverage,FindMaxMin and FindPercentage */ /* functions to fix a potential bug. */ /* 1.0 30.11.1995 First official release */ /****************************************************************************/ #include #include #include #include "NewTime.h" #include "StatInfo.h" #define MissingValue 9999999 /* Marker for missing data value */ #define OK 0 /* Return value for successfull operation */ #define FAIL 1 /* Return value for failed operation */ #define FileError 1 /* Error code for failed file operation */ #define OutOfMemory 2 /* Error code for memory allocation failure */ #define FileFormatError 3 /* Error code for wrong file format type */ #define StationNotFound 4 /* Station not found in the data file */ /*--------------------------------------------------------------------------*/ /* Here is a structure which contains data for one day. The strategy for */ /* reading data from a file or from stdin is that one first reads the data */ /* into a chained list of OneDay structures. After all the data has been */ /* read the data is moved into proper place in the Station_struct. This */ /* algorithm might sound complex but is necessary since one doesn't know */ /* a priori the amount of data to be read (e.g. if it comes from stdin). */ /* Therefore the memory allocation must be done in small pieces. */ /*--------------------------------------------------------------------------*/ struct OneDay { /* Data for one observatory for one day */ long *X,*Y,*Z; /* Pointers to data areas of the first */ /* three components (usually X, Y and Z). */ long *H,*D,*F; /* Pointers to data areas of three more */ /* components (e.g. H, D and F). */ long *T; /* Pointer to Temperature data (T). */ long *A,*B; /* Pointers to extra A and B components */ /* that are used e.g. in Dump-format files. */ long XCount,YCount,ZCount; /* Number of datapoints in data blocks */ /* of X,Y,Z,H,D and F components. */ long TCount; /* Number of datapoints in data blocks */ /* of T,A and B components. */ /* Number of A and B components are assumed */ struct OneDay *NextDay; /* Pointer to the next one day block */ }; typedef struct OneDay OneDay_struct; typedef OneDay_struct *OneDayPtr; /*--------------------------------------------------------------------------*/ /* The data structure for one recording station. */ /* The unit of time is second. The units used for magnetic field components */ /* are defined by the user. Routines defined here do not assume anything */ /* about the magnitude of the field values. Typically the units are defined */ /* in the routines which read data in (e.g. IAGA.h, GADF.h, WDC.h). One */ /* must then be careful to convert the units when e.g. reading and writing */ /* data in different formats (e.g. conversion Dump -> IAGA). */ /* Note, however, that long's are used to represent magnetic field values. */ /* Typical unit for field value is 0.1 nT (IAGA.h) or 0.01 nT (Dump.h). */ /* Unit for temperature is similarly undefined. Units for geographical */ /* coordinates (latitude and longitude) is 0.01 degrees. */ /*--------------------------------------------------------------------------*/ struct Station { char StationID[4]; /* Three letter station ID */ long Latitude; /* Geographic latitude of the station (0.01 degs) */ long Longitude; /* Geographic longitude of the station (0.01 degs) */ Time_sec StartTime; /* Time of the first data point */ Time_sec EndTime; /* Time of the last data point + TimeStep */ char Components[9]; /* Components: e.g. "XYZ" or "HDZ" */ long *X,*Y,*Z; /* Pointers to data areas of the first three */ /* components (usually X, Y and Z). */ long *H,*D,*F; /* Pointers to data areas of three more */ /* components (e.g. H, D and F). */ long *T; /* Pointer to Temperature data (T). */ long *A,*B; /* Pointers to extra A and B components that are */ /* used e.g. in Dump-format files. */ long TimeStep; /* Time difference between successive data points */ /* of X,Y,Z,H,D and F components. */ long TempStep; /* Time separation between successive data points */ /* of T,A and B components. */ struct Station *Next; /* Pointer to the next station in the list */ OneDayPtr FirstDay; /* Pointer to the first one day structure. */ /* This pointer is used only during data loading, */ /* after that it is set to nil value. */ long Xbase; /* Extra space for storing baseline, average etc. */ long Ybase; /* Extra space for storing baseline, average etc. */ long Zbase; /* Extra space for storing baseline, average etc. */ long Hbase; /* Extra space for storing baseline, average etc. */ long Dbase; /* Extra space for storing baseline, average etc. */ long Fbase; /* Extra space for storing baseline, average etc. */ }; typedef struct Station Station_struct; typedef Station_struct *StationPtr; /*--------------------------------------------------------------------------*/ /* The data structure for a network of recording stations */ /*--------------------------------------------------------------------------*/ struct Network { char Name[20]; /* Name of the network */ long StationCount; /* Number of stations in the network */ StationPtr StationList; /* Pointer to the first station */ }; typedef struct Network Network_struct; typedef Network_struct *NetworkPtr; /*--------------------------------------------------------------------------*/ /* Function prototypes */ /*--------------------------------------------------------------------------*/ void InitNetwork(NetworkPtr Network); void FreeNetwork(NetworkPtr Network); StationPtr FindStation(NetworkPtr Network, char *StationName); long StationInList(char *Station, char *StationList); StationInfoPtr FindStationInfo(char *StationID); static void FreeOneDayBlock(OneDayPtr Day); void FillWithMissingData(long *DataPtr,long Count); void CopyDataBlock(long *From, long *To, long Count); static OneDayPtr NewOneDayBlock(long DataStep,long TempStep); long OneDayFull(StationPtr Station,char Component); long AddOneDay(StationPtr Station); StationPtr AddNewStation(NetworkPtr Network,long DataStep,long TempStep); long NewStation(NetworkPtr Network,StationInfoPtr p,Time_sec StartTime, Time_sec EndTime, Time_sec TimeStep, Time_sec TempStep); void RemoveStation(NetworkPtr Network,char *StationName); long CombineData(StationPtr Station); long *GetDataPtr(StationPtr Station,Time_sec Secs,char Comp); long GetDataValue(StationPtr Station,Time_sec Secs,char Comp); void SetDataValue(StationPtr s,Time_sec Secs,char Comp,long Value); long ComputeAverage(StationPtr Station,char Comp,Time_sec Time,long Length); long ComputeAverage_Conditional(StationPtr Station,char Comp,Time_sec Time, long Length, long Percentage); void FindMaxMin(StationPtr Station,char Comp,Time_sec Time,long Length, long *MaxValue, long *MinValue); long FindMaxAplitude(StationPtr Station,Time_sec StartTime,long Length); long FindPercentage(StationPtr Station,char Comp,Time_sec Time,long Length); long RoundFloat(double x); long FindGap(StationPtr Station, char Component, Time_sec T, Time_sec *GapStart, Time_sec *GapEnd); void InterpolateData(StationPtr Station,Time_sec T1, Time_sec T2, char Comp); void CopyMagnData(StationPtr FromStation, StationPtr ToStation, char Comp, Time_sec StartTime, Time_sec EndTime, long BaseShift); void CopyAllMagnData(StationPtr FromStation, StationPtr ToStation, Time_sec StartTime, Time_sec EndTime, long BaseShift); void ExtendDataArea(StationPtr s, Time_sec StartTime, Time_sec EndTime); void ChangeUnitsUP(NetworkPtr Network); void ChangeUnitsDOWN(NetworkPtr Network); void ScaleData(StationPtr Station, char Comp, double Coeff); /*--------------------------------------------------------------------------*/ /* Initialize a Network structure */ /*--------------------------------------------------------------------------*/ void InitNetwork(NetworkPtr Network) { Network->Name[0] = '\0'; Network->StationCount = 0; Network->StationList = NULL; } /*--------------------------------------------------------------------------*/ /* Release the space allocated for data and stations in a Network structure */ /*--------------------------------------------------------------------------*/ void FreeNetwork(NetworkPtr Network) { StationPtr s,r; OneDayPtr p,q; for (s = Network->StationList; s != NULL; s = r) { if ((s->X) != NULL) free(s->X); /* Release data areas */ if ((s->Y) != NULL) free(s->Y); if ((s->Z) != NULL) free(s->Z); if ((s->H) != NULL) free(s->H); if ((s->D) != NULL) free(s->D); if ((s->F) != NULL) free(s->F); if ((s->T) != NULL) free(s->T); if ((s->A) != NULL) free(s->A); if ((s->B) != NULL) free(s->B); /* Release OneDayBlocks */ for (p = s->FirstDay; p != NULL; p = q) { q = p->NextDay; FreeOneDayBlock(p); } r = s->Next; free(s); /* Release the space of station structure */ } InitNetwork(Network); } /*--------------------------------------------------------------------------*/ /* Try to find a station from the Network. If found return the pointer to */ /* the station structure, otherwise return NULL. */ /*--------------------------------------------------------------------------*/ StationPtr FindStation(NetworkPtr Network, char *StationName) { StationPtr s; s = Network->StationList; while ((s != NULL) && (strncmp(s->StationID,StationName,3))) s = s->Next; return (s); } /*--------------------------------------------------------------------------*/ /* Check if the three character station ID 'Station' is in the list of */ /* stations. The StationList is a string with 3-letter station ID's */ /* separated by a space (e.g. 'SOR MAS HAN NUR'). If the StationList is */ /* NULL or zero length then true value ( = 1) is returned. */ /*--------------------------------------------------------------------------*/ long StationInList(char *Station, char *StationList) { long found; char *q = StationList; if ((q == NULL) || (q[0] == '\0')) return 1; do { /* Go through all 3-letter ID's */ found = (strncmp(Station,q,3) == 0); q += 3; } while ((*q++ == ' ') && (!found)); return (found); } /*--------------------------------------------------------------------------*/ /* Search the StationInfoTable defined in the header file 'StatInfo.h' and */ /* find a pointer to the StationInfoStruct corresponding to the given */ /* station. If station is not found return a pointer to the last element */ /* of StationInfoTable which contains some default values. */ /*--------------------------------------------------------------------------*/ StationInfoPtr FindStationInfo(char *StationID) { long i = 0; StationInfoPtr p = StationInfoTable; /* Compare the StationID with the names defined in the StationInfoTable */ while ((++i < IMAGEStationCount) && (strncmp(p->StationID,StationID,3))) p++; return (p); } /*--------------------------------------------------------------------------*/ /* Release the area allocated for a OneDayBlock and for the accompanying */ /* data. */ /*--------------------------------------------------------------------------*/ static void FreeOneDayBlock(OneDayPtr Day) { if ((Day->X) != NULL) free(Day->X); if ((Day->Y) != NULL) free(Day->Y); if ((Day->Z) != NULL) free(Day->Z); if ((Day->H) != NULL) free(Day->H); if ((Day->D) != NULL) free(Day->D); if ((Day->F) != NULL) free(Day->F); if ((Day->T) != NULL) free(Day->T); if ((Day->A) != NULL) free(Day->A); if ((Day->B) != NULL) free(Day->B); free(Day); } /*--------------------------------------------------------------------------*/ /* Fill given data area with missing value markers. */ /*--------------------------------------------------------------------------*/ void FillWithMissingData(long *DataPtr,long Count) { long i; long *p = DataPtr; for (i=0;iX = (long *) malloc(BlockSizeData); p->Y = (long *) malloc(BlockSizeData); p->Z = (long *) malloc(BlockSizeData); p->H = (long *) malloc(BlockSizeData); p->D = (long *) malloc(BlockSizeData); p->F = (long *) malloc(BlockSizeData); p->T = (long *) malloc(BlockSizeTemp); p->A = (long *) malloc(BlockSizeTemp); p->B = (long *) malloc(BlockSizeTemp); p->XCount = 0; p->YCount = 0; p->ZCount = 0; p->TCount = 0; p->NextDay = NULL; if (((p->X) == NULL) || ((p->Y) == NULL) || ((p->Z) == NULL) || ((p->H) == NULL) || ((p->D) == NULL) || ((p->F) == NULL) || ((p->T) == NULL) || ((p->A) == NULL) || ((p->B) == NULL)) { FreeOneDayBlock(p); return (NULL); } FillWithMissingData(p->X,86400/DataStep); FillWithMissingData(p->Y,86400/DataStep); FillWithMissingData(p->Z,86400/DataStep); FillWithMissingData(p->H,86400/DataStep); FillWithMissingData(p->D,86400/DataStep); FillWithMissingData(p->F,86400/DataStep); FillWithMissingData(p->T,86400/TempStep); FillWithMissingData(p->A,86400/TempStep); FillWithMissingData(p->B,86400/TempStep); return(p); } /*--------------------------------------------------------------------------*/ /* Check if the last OneDayBlock is already full of data for given */ /* component. */ /*--------------------------------------------------------------------------*/ long OneDayFull(StationPtr Station,char Component) { OneDayPtr p; long Count; /* Find the last OneDayBlock */ for (p = Station->FirstDay ; p->NextDay != NULL ; p = p->NextDay); switch (Component) { case 'H' : case 'X' : Count = p->XCount; break; case 'D' : case 'Y' : Count = p->YCount; break; case 'Z' : Count = p->ZCount; break; } return (Count == ((long) 86400)/(Station->TimeStep)); } /*--------------------------------------------------------------------------*/ /* Add a new OneDayBlock into the Station. If memory allocation isn't */ /* successfull return 1, otherwise return 0. */ /*--------------------------------------------------------------------------*/ long AddOneDay(StationPtr Station) { OneDayPtr p; /* Find the last OneDayBlock */ for (p = Station->FirstDay ; p->NextDay != NULL ; p = p->NextDay); /* Create a new OneDayBlock and set the pointers */ p->NextDay = NewOneDayBlock(Station->TimeStep,Station->TempStep); if ((p->NextDay) == NULL) return 1; else return 0; } /*--------------------------------------------------------------------------*/ /* Create a new Station_struct and create also the first OneDayBlock. */ /* Note however that most of the fields in the Station_struct are left */ /* undefined. If memory allocation fails then return a NULL pointer. */ /*--------------------------------------------------------------------------*/ StationPtr AddNewStation(NetworkPtr Network,long DataStep,long TempStep) { StationPtr s,q; s = (StationPtr) malloc(sizeof(Station_struct)); if (s == NULL) return(NULL); s->TimeStep = DataStep; s->TempStep = TempStep; strcpy(s->Components,"???"); s->X = NULL; s->Y = NULL; s->Z = NULL; s->H = NULL; s->D = NULL; s->F = NULL; s->T = NULL; s->A = NULL; s->B = NULL; s->Next = NULL; s->FirstDay = NewOneDayBlock(DataStep,TempStep); /* Add the station to the Network */ if ((Network->StationList) == NULL) Network->StationList = s; else { q = Network->StationList; while ((q->Next) != NULL) q = q->Next; q->Next = s; } (Network->StationCount)++; if ((s->FirstDay) == NULL) return (NULL); else return (s); } /*--------------------------------------------------------------------------*/ /* Gather the data from the list of OneDayBlocks into a continuous data */ /* area and release the areas allocated to OneDayBlocks. If no memory can */ /* be allocated OutOfMemory error will be returned, otherwise OK. */ /*--------------------------------------------------------------------------*/ long CombineData(StationPtr s) { long CountData; long CountTemp; long BlockSizeTemp; long BlockSizeData; long *pX,*pY,*pZ,*pH,*pD,*pF,*pT,*pA,*pB; OneDayPtr r,q; /* Find the number of datapoints and allocate space for them */ CountData = (s->EndTime - s->StartTime)/(s->TimeStep); CountTemp = (s->EndTime - s->StartTime)/(s->TempStep); BlockSizeData = sizeof(long)*CountData; BlockSizeTemp = sizeof(long)*CountTemp; s->X = (long *) malloc(BlockSizeData); s->Y = (long *) malloc(BlockSizeData); s->Z = (long *) malloc(BlockSizeData); s->H = (long *) malloc(BlockSizeData); s->D = (long *) malloc(BlockSizeData); s->F = (long *) malloc(BlockSizeData); s->T = (long *) malloc(BlockSizeTemp); s->A = (long *) malloc(BlockSizeTemp); s->B = (long *) malloc(BlockSizeTemp); if ((s->X == NULL) || (s->Y == NULL) || (s->Z == NULL) || (s->H == NULL) || (s->D == NULL) || (s->F == NULL) || (s->T == NULL) || (s->A == NULL) || (s->B == NULL)) return (OutOfMemory); /* Fill data blocks with missing data */ FillWithMissingData(s->X,CountData); FillWithMissingData(s->Y,CountData); FillWithMissingData(s->Z,CountData); FillWithMissingData(s->H,CountData); FillWithMissingData(s->D,CountData); FillWithMissingData(s->F,CountData); FillWithMissingData(s->T,CountTemp); FillWithMissingData(s->A,CountTemp); FillWithMissingData(s->B,CountTemp); /* Copy data from the OneDayBlocks */ pX = s->X; pY = s->Y; pZ = s->Z; pH = s->H; pD = s->D; pF = s->F; pT = s->T; pA = s->A; pB = s->B; for (r = s->FirstDay;r != NULL;r = q) { memcpy(pX,r->X,sizeof(long)*(r->XCount)); pX += r->XCount; memcpy(pY,r->Y,sizeof(long)*(r->YCount)); pY += r->YCount; memcpy(pZ,r->Z,sizeof(long)*(r->ZCount)); pZ += r->ZCount; memcpy(pH,r->H,sizeof(long)*(r->ZCount)); pH += r->ZCount; memcpy(pD,r->D,sizeof(long)*(r->ZCount)); pD += r->ZCount; memcpy(pF,r->F,sizeof(long)*(r->ZCount)); pF += r->ZCount; memcpy(pT,r->T,sizeof(long)*(r->TCount)); pT += r->TCount; memcpy(pA,r->A,sizeof(long)*(r->TCount)); pA += r->TCount; memcpy(pB,r->B,sizeof(long)*(r->TCount)); pB += r->TCount; q = r->NextDay; FreeOneDayBlock(r); } s->FirstDay = NULL; return OK; } /*--------------------------------------------------------------------------*/ /* Add station to the Network structure and fill the parameters in the */ /* Station structure for this particular station. Also fill the data */ /* values with missing data markers */ /*--------------------------------------------------------------------------*/ long NewStation(NetworkPtr Network,StationInfoPtr p,Time_sec StartTime, Time_sec EndTime, Time_sec TimeStep, Time_sec TempStep) { StationPtr s,q; long j,CountData,CountTemp; long BlockSizeTemp,BlockSizeData; long *pX,*pY,*pZ,*pH,*pD,*pF,*pT,*pA,*pB; q = Network->StationList; if ((q != NULL) && (StartTime == 0)) StartTime = q->StartTime; if ((q != NULL) && (EndTime == 0)) EndTime = q->EndTime; if ((q != NULL) && (TimeStep == 0)) TimeStep = q->TimeStep; if ((q != NULL) && (TempStep == 0)) TempStep = q->TempStep; /* Fill some parameters */ s = (StationPtr) malloc(sizeof(Station_struct)); strcpy(s->StationID,p->StationID); s->Latitude = p->Latitude; s->Longitude = p->Longitude; s->StartTime = StartTime; s->EndTime = EndTime; s->TimeStep = TimeStep; s->TempStep = TempStep; strcpy(s->Components,"XYZ"); s->Next = NULL; s->FirstDay = NULL; /* Add the station to the Network */ if (q == NULL) Network->StationList = s; else { while ((q->Next) != NULL) q = q->Next; q->Next = s; } (Network->StationCount)++; /* Find the number of datapoints and allocate space for them */ CountData = (EndTime - StartTime)/TimeStep; CountTemp = (EndTime - StartTime)/TempStep; BlockSizeData = sizeof(long)*CountData; BlockSizeTemp = sizeof(long)*CountTemp; s->X = (long *) malloc(BlockSizeData); s->Y = (long *) malloc(BlockSizeData); s->Z = (long *) malloc(BlockSizeData); s->H = (long *) malloc(BlockSizeData); s->D = (long *) malloc(BlockSizeData); s->F = (long *) malloc(BlockSizeData); s->T = (long *) malloc(BlockSizeTemp); s->A = (long *) malloc(BlockSizeTemp); s->B = (long *) malloc(BlockSizeTemp); if ((s->X == NULL) || (s->Y == NULL) || (s->Z == NULL) || (s->H == NULL) || (s->D == NULL) || (s->F == NULL) || (s->T == NULL) || (s->A == NULL) || (s->B == NULL)) return (OutOfMemory); s->Xbase = MissingValue; s->Ybase = MissingValue; s->Zbase = MissingValue; s->Hbase = MissingValue; s->Dbase = MissingValue; s->Fbase = MissingValue; /* Fill the data fields with missing values */ pX = s->X; pY = s->Y; pZ = s->Z; pH = s->H; pD = s->D; pF = s->F; pT = s->T; pA = s->A; pB = s->B; for (j=0;jStationList; while ((s != NULL) && (strncmp(s->StationID,StationName,3))) { p = s; s = s->Next; } /* Remove the station from the linked list and free the data space */ if (s != NULL) { if (s == Network->StationList) /* First station in the list */ Network->StationList = s->Next; else p->Next = s->Next; Network->StationCount--; /* Decrease the number of stations */ if ((s->X) != NULL) free(s->X); /* Release data areas */ if ((s->Y) != NULL) free(s->Y); if ((s->Z) != NULL) free(s->Z); if ((s->H) != NULL) free(s->H); /* Release data areas */ if ((s->D) != NULL) free(s->D); if ((s->F) != NULL) free(s->F); if ((s->T) != NULL) free(s->T); if ((s->A) != NULL) free(s->A); if ((s->B) != NULL) free(s->B); free(s); /* Release the space of station structure */ } } /*--------------------------------------------------------------------------*/ /* Get a pointer to a data value for specified station, time and component. */ /* If there is no data for the specified time a NULL pointer is returned. */ /*--------------------------------------------------------------------------*/ long *GetDataPtr(StationPtr s,Time_sec Secs,char Comp) { long *DataPtr; long count; /* Check the times */ if ((Secs < s->StartTime) || (Secs >= s->EndTime) || (((Secs - s->StartTime) % s->TimeStep) != 0)) return (NULL); /* Get the component */ switch (Comp) { case 'X' : DataPtr = s->X; break; case 'Y' : DataPtr = s->Y; break; case 'Z' : DataPtr = s->Z; break; case 'H' : DataPtr = s->H; break; case 'D' : DataPtr = s->D; break; case 'F' : DataPtr = s->F; break; case 'T' : DataPtr = s->T; break; case 'A' : DataPtr = s->A; break; case 'B' : DataPtr = s->B; break; default : return (NULL); } if (DataPtr == NULL) return (NULL); /* Finally get the pointer */ if ((Comp == 'T') || (Comp == 'A') || (Comp == 'B')) count = (Secs - s->StartTime)/(s->TempStep); else count = (Secs - s->StartTime)/(s->TimeStep); return (DataPtr+count); } /*--------------------------------------------------------------------------*/ /* Get the data value for specified station, time and component. If there */ /* is no data for the specified time a missing value is returned. */ /* This routine is not very fast because of various checkings. If one wants */ /* to retrieve consecutive data points it is much faster to get a pointer */ /* to the first data point with GetDataPtr (see above) and then use that */ /* pointer to retrieve the consecutive data values. */ /*--------------------------------------------------------------------------*/ long GetDataValue(StationPtr s,Time_sec Secs,char Comp) { long *DataPtr; DataPtr = GetDataPtr(s,Secs,Comp); if (DataPtr != NULL) return(*DataPtr); else return(MissingValue); } /*--------------------------------------------------------------------------*/ /* Set the data value for specified station, time and component. If the */ /* station data does not include the specified time point then nothing is */ /* done. */ /* This routine is not very fast because of various checkings. If one wants */ /* to set consecutive data points it is much faster to get a pointer to the */ /* first data point with GetDataPtr (see above) and then use that pointer */ /* to set consecutive data values. */ /*--------------------------------------------------------------------------*/ void SetDataValue(StationPtr s,Time_sec Secs,char Comp,long Value) { long *DataPtr; DataPtr = GetDataPtr(s,Secs,Comp); if (DataPtr != NULL) *DataPtr = Value; } /*--------------------------------------------------------------------------*/ /* Round a floating point number into a long. Coudn't find it in any of */ /* the standard libraries !?. */ /*--------------------------------------------------------------------------*/ long RoundFloat(double x) { if (x < 0) return ((long)(x-0.5)); else return ((long)(x+0.5)); } /*--------------------------------------------------------------------------*/ /* Compute an average over speficied time period for specified station, */ /* component and time. */ /* An average is computed even if a single valid data point exist in the */ /* given interval. */ /*--------------------------------------------------------------------------*/ long ComputeAverage(StationPtr Station,char Comp,Time_sec Time,long Length) { double Sum = 0; /* we must use double as long might overflow */ long Good = 0; long i,Count,Value; if (Length == 0) Count = 1; /* Average is same as a single data value */ else Count = Length/Station->TimeStep; for (i=0;iTimeStep; } if (Good == 0) return (MissingValue); else return RoundFloat(Sum/Good); } /*--------------------------------------------------------------------------*/ /* Compute an average over speficied time period for specified station, */ /* component and time. */ /* This is same as above but with the condition that a minimum percentage */ /* of valid data points must exist in the given interval before an average */ /* can be computed. */ /*--------------------------------------------------------------------------*/ long ComputeAverage_Conditional(StationPtr Station,char Comp,Time_sec Time, long Length, long Percentage) { double Sum = 0; /* we must use double as long might overflow */ long Good = 0; long i,Count,Value; if (Length == 0) Count = 1; /* Average is same as a single data value */ else Count = Length/Station->TimeStep; for (i=0;iTimeStep; } if (Good == 0) return (MissingValue); if (Good < (Count*Percentage)/100) return (MissingValue); return RoundFloat(Sum/Good); } /*--------------------------------------------------------------------------*/ /* Find the maximum and minimum field values over speficied time period */ /* station, component and time. */ /*--------------------------------------------------------------------------*/ void FindMaxMin(StationPtr Station,char Comp,Time_sec Time,long Length, long *MaxValue, long *MinValue) { long i,Count,Value; long Good = 0; long Max; long Min; if (Length == 0) Count = 1; /* Max and Min are same */ else Count = Length/Station->TimeStep; for (i=0;i Max) Max = Value; else if (Value < Min) Min = Value; } Time += Station->TimeStep; } if (Good == 0) { *MaxValue = MissingValue; *MinValue = MissingValue; } else { *MaxValue = Max; *MinValue = Min; } } /*--------------------------------------------------------------------------*/ /* Find the percentage of non-missing data points in the given time period,*/ /* station and component. */ /*--------------------------------------------------------------------------*/ long FindPercentage(StationPtr Station,char Comp,Time_sec Time,long Length) { long i,Count,Value; long Good = 0; if (Length == 0) Count = 1; /* Average is same as a single data value */ else Count = Length/Station->TimeStep; for (i=0;iTimeStep; } return (100*Good)/Count; } /*--------------------------------------------------------------------------*/ /* Find the next data gap for given station and component starting from */ /* time T. If no gap is found 0 is returned and GapEnd is set to */ /* Station->EndTime. If gap is found GapStart is set to the time of first */ /* missing data point and GapEnd is set to the time of first existing data */ /* point after missing value (=last missing data value+Station->TimeStep). */ /*--------------------------------------------------------------------------*/ long FindGap(StationPtr Station, char Component, Time_sec T, Time_sec *GapStart, Time_sec *GapEnd) { long *DataPtr; DataPtr = GetDataPtr(Station,T,Component); /*** Find the first missing data point ***/ while ((T < Station->EndTime) && (*DataPtr != MissingValue)) { DataPtr++; T += Station->TimeStep; } if (T == Station->EndTime) { /* No missing value found */ *GapEnd = Station->EndTime; return 0; } /*** Find the last missing data point ***/ *GapStart = T; while ((T < Station->EndTime) && (*DataPtr == MissingValue)) { DataPtr++; T += Station->TimeStep; } *GapEnd = T; return 1; } /*--------------------------------------------------------------------------*/ /* Interpolate data points linearly from time 'T1' to time 'T2' based on */ /* field values at 'T1-TimeStep' and 'T2+TimeStep'. If one of these is */ /* a missing value then extrapolation is used. If both are missing values */ /* then also the interpolated data points will be set as missing. */ /*--------------------------------------------------------------------------*/ void InterpolateData(StationPtr Station,Time_sec T1, Time_sec T2, char Comp) { long X1,X2,Value; Time_sec T,Step; double coeff; if ((Comp == 'T') || (Comp == 'A')) Step = Station->TempStep; else Step = Station->TimeStep; X1 = GetDataValue(Station,T1 - Step,Comp); X2 = GetDataValue(Station,T2 + Step,Comp); /* If both end points are missing then interpolate with missing value */ if ((X1 == MissingValue) && (X2 == MissingValue)) { for (T=T1; T<=T2; T += Step) SetDataValue(Station,T,Comp,MissingValue); return; } /* One of the end points is not missing, therefore either */ /* interpolate or extrapolate. */ if (X1 == MissingValue) X1 = X2; if (X2 == MissingValue) X2 = X1; coeff = ((double)(X2-X1))/((double)(T2-T1+2*Step)); for (T=T1; T<=T2; T += Step) { Value = RoundFloat(X1+coeff*(T-T1+Step)); SetDataValue(Station,T,Comp,Value); } } /*--------------------------------------------------------------------------*/ /* Copy data from one station structure to another. Also adjust for the */ /* possible baseline shift. */ /*--------------------------------------------------------------------------*/ void CopyMagnData(StationPtr FromStation, StationPtr ToStation,char Comp, Time_sec StartTime, Time_sec EndTime, long BaseShift) { Time_sec Time,Step; long *DataPtr; long Value; if ((Comp == 'T') || (Comp == 'A') || (Comp == 'B')) Step = FromStation->TempStep; else Step = FromStation->TimeStep; Time = StartTime; while (Time < EndTime) { DataPtr = GetDataPtr(ToStation,Time,Comp); if (DataPtr != NULL) { Value = GetDataValue(FromStation,Time,Comp); if (Value != MissingValue) *DataPtr = Value+BaseShift; else *DataPtr = MissingValue; } Time += Step; } } void CopyAllMagnData(StationPtr FromStation, StationPtr ToStation, Time_sec StartTime, Time_sec EndTime, long BaseShift) { char Comps[10] = "XYZHDFTAB"; long i; for (i=0;i<9;i++) CopyMagnData(FromStation,ToStation, Comps[i], StartTime, EndTime, BaseShift); } /*--------------------------------------------------------------------------*/ /* Enlarge the data area of station data to given time limits. The added */ /* space will be filled with missing data. */ /*--------------------------------------------------------------------------*/ void ExtendDataArea(StationPtr s, Time_sec StartTime, Time_sec EndTime) { Time_sec T0 = (StartTime < s->StartTime) ? StartTime : s->StartTime; Time_sec T1 = (EndTime > s->EndTime) ? EndTime : s->EndTime; long CountData; long CountTemp; long BlockSizeData; long BlockSizeTemp; long OffsetData; long OffsetTemp; long *X,*Y,*Z,*H,*D,*F,*T,*A,*B; /* Find the number of datapoints and allocate space for them */ CountData = (T1 - T0)/s->TimeStep; CountTemp = (T1 - T0)/s->TempStep; BlockSizeData = sizeof(long)*CountData; BlockSizeTemp = sizeof(long)*CountTemp; X = (long *) malloc(BlockSizeData); Y = (long *) malloc(BlockSizeData); Z = (long *) malloc(BlockSizeData); H = (long *) malloc(BlockSizeData); D = (long *) malloc(BlockSizeData); F = (long *) malloc(BlockSizeData); T = (long *) malloc(BlockSizeTemp); A = (long *) malloc(BlockSizeTemp); B = (long *) malloc(BlockSizeTemp); /* Fill the data fields with missing values */ FillWithMissingData(X,CountData); FillWithMissingData(Y,CountData); FillWithMissingData(Z,CountData); FillWithMissingData(H,CountData); FillWithMissingData(D,CountData); FillWithMissingData(F,CountData); FillWithMissingData(T,CountTemp); FillWithMissingData(A,CountTemp); FillWithMissingData(B,CountTemp); /* Copy data from old data areas */ CountData = (s->EndTime - s->StartTime)/s->TimeStep; CountTemp = (s->EndTime - s->StartTime)/s->TempStep; OffsetData = (s->StartTime - T0)/s->TimeStep; OffsetTemp = (s->StartTime - T0)/s->TempStep; CopyDataBlock(s->X,X+OffsetData,CountData); CopyDataBlock(s->Y,Y+OffsetData,CountData); CopyDataBlock(s->Z,Z+OffsetData,CountData); CopyDataBlock(s->H,H+OffsetData,CountData); CopyDataBlock(s->D,D+OffsetData,CountData); CopyDataBlock(s->F,F+OffsetData,CountData); CopyDataBlock(s->T,T+OffsetTemp,CountTemp); CopyDataBlock(s->A,A+OffsetTemp,CountTemp); CopyDataBlock(s->B,B+OffsetTemp,CountTemp); /* Release previous data areas */ if ((s->X) != NULL) free(s->X); if ((s->Y) != NULL) free(s->Y); if ((s->Z) != NULL) free(s->Z); if ((s->H) != NULL) free(s->H); if ((s->D) != NULL) free(s->D); if ((s->F) != NULL) free(s->F); if ((s->T) != NULL) free(s->T); if ((s->A) != NULL) free(s->A); if ((s->B) != NULL) free(s->B); /* Update the pointers */ s->X = X; s->Y = Y; s->Z = Z; s->H = H; s->D = D; s->F = F; s->T = T; s->A = A; s->B = B; s->StartTime = T0; s->EndTime = T1; } /*--------------------------------------------------------------------------*/ /* Multiply all field values by 10. This corresponds to changing the field */ /* units from 1 nT to 0.1 nT or from 0.1 nT to 0.01 nT. */ /*--------------------------------------------------------------------------*/ static void Multiply10(long *DataPtr, long count) { long i; for (i=0;iStationList; while (s != NULL) { /* Go through all components */ count = (s->EndTime - s->StartTime)/s->TimeStep; Multiply10(s->X,count); Multiply10(s->Y,count); Multiply10(s->Z,count); Multiply10(s->H,count); Multiply10(s->D,count); Multiply10(s->F,count); count = (s->EndTime - s->StartTime)/s->TempStep; Multiply10(s->T,count); s = s->Next; } } /*--------------------------------------------------------------------------*/ /* Change the magnetic field units from 0.1 nT to 1 nT or from 0.01 nT to */ /* 0.1 nT. This is done by dividing all field values by 10. */ /* Also change temperature units from 0.01 C to 0.1 Centigrades. */ /*--------------------------------------------------------------------------*/ void ChangeUnitsDOWN(NetworkPtr Network) { StationPtr s; long count; /* Go through all the stations */ s = Network->StationList; while (s != NULL) { /* Go through all components */ count = (s->EndTime - s->StartTime)/s->TimeStep; Divide10(s->X,count); Divide10(s->Y,count); Divide10(s->Z,count); Divide10(s->H,count); Divide10(s->D,count); Divide10(s->F,count); count = (s->EndTime - s->StartTime)/s->TempStep; Divide10(s->T,count); s = s->Next; } } /*--------------------------------------------------------------------------*/ /* Scale data by given coefficient. */ /*--------------------------------------------------------------------------*/ void ScaleData(StationPtr Station, char Comp, double Coeff) { Time_sec T; long Value; for (T=Station->StartTime; TEndTime; T += Station->TimeStep) { Value = GetDataValue(Station,T,Comp); if (Value != MissingValue) { Value = RoundFloat(Coeff*Value); SetDataValue(Station,T,Comp,Value); } } }