/****************************************************************************/
/*                                                                          */
/*                                WDChour.h                                 */
/*                                                                          */
/****************************************************************************/
/****************************************************************************/
/* This is a C language header file that defines the following routines     */
/* for reading one hour WDC-format data files.                              */
/*                                                                          */
/*  ----------------  Read data from a one hour WDC-format file ----------  */
/*                                                                          */
/* long ReadWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,   */
/*               char *EndTimeStr,char *StationList)                        */
/*      FileName        Name of the WDC-file containing the data.           */
/*                      stdin in not allowed.                               */
/*      Network         Structure into which the data is read.              */
/*                      Network_struct is defined in file MagnData.h.       */
/*      StartTimeStr    String defining the first time to be read from      */
/*                      the file (String delimited with '\0').              */
/*                      String format : YYMMDDHH. If HH is missing then     */
/*                      00 is assumed. If StartTimeStr is nil or zero       */
/*                      length then start of file is assumed.               */
/*      EndTimeStr      String defining the time not to be read anymore     */
/*                      Format as before YYMMDDHH. If nil or zero length    */
/*                      then data is read until end of file.                */
/*      StationList     List of stations included in the data. The stations */
/*                      must be specified with three letter codes (e.g.     */
/*                      SOR) and must be separeted with a space in the list */
/*                      (e.g. 'SOR MAS KEV'). If nil or zero length then    */
/*                      all stations in the file are be included.           */
/*                                                                          */
/*  --------------  Write data into a one hour WDC-format file  ----------  */
/*                                                                          */
/* long WriteWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,  */
/*                char *EndTimeStr,char *StationList)                       */
/*      FileName        Name of the WDC file. IF NULL or zero length then   */
/*                      stdout is used.                                     */
/*      Network         A pointer to the Network structure whose data is to */
/*                      be written into a file.                             */
/*      StartTime       Start date string YYMMDDHHMM. If NULL then same as  */
/*                      the time of the first record of the first station   */
/*                      in the Network data structure.                      */
/*      EndTime         End date string YYMMDDHHMM. This record is NOT      */
/*                      included anymore. If NULL then data till the end of */
/*                      available data in Network structure will be written */
/*      StationList     String containing the 3-letter ID's of the stations */
/*                      to be included. If NULL or 0-length then all the    */
/*                      stations will be written. The stations must be      */
/*                      separated by a space in the list (e.g. 'SOR MUO').  */
/****************************************************************************/
/****************************************************************************/
/*                            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.04    12.02.2020                          */
/****************************************************************************/
/****************************************************************************/
/*  Version history:                                                        */
/*                                                                          */
/*  1.04 12.02.2020 Date strings must be YYYYMMDD (four digits for year)    */
/*                  Replaced StrToSecs -> StrToSecsC and                    */
/*                  SecsToStr -> SecsToStrC.                                */
/*  1.03 03.11.2014  Added HDZ <-> XYZ conversions                          */
/*  1.02 10.11.2000  Total rewrite of the ReadWDChour routine               */
/*  1.01 08.11.2000  Added WriteWDChour routine                             */
/*  1.0  24.10.2000  First release                                          */
/****************************************************************************/


#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>

#define WDCstep 3600        /* One-hour values. No other possibilities !    */
#define MissingWDC 9999     /* Marker for missing data value                */
#define WDCLineLength 122   /* Length of one data line (including linefeed) */

/*--------------------------------------------------------------------------*/
/* Define the offsets of various fields from the start of the WDC block.    */
/* See the definition of WDC record for the explanation of various fields.  */
/*--------------------------------------------------------------------------*/

#define WDC_Station      0
#define WDC_YYMM         3
#define WDC_Component    7
#define WDC_Day          8
#define WDC_Century     14
#define WDC_Base        16
#define WDC_Data        20
#define WDC_DailyMean  116

/*--------------------------------------------------------------------------*/
/* Function prototypes                                                      */
/*--------------------------------------------------------------------------*/

static long GetFileSize(char *FileName);
static Time_sec GetWDCTime(char *Buffer);
static void SetStationParamsWDC(StationPtr Station,char *Buff);
static long GetWDCFieldValue(char *p);
static void GetWDCData(StationPtr Station, Time_sec Time, char *Buff);
long ReadWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,
              char *EndTimeStr,char *StationList, long *Century);

static void WriteWDCBlock(FILE *DataFile,StationPtr Station,Time_sec Time,
                   char Component);
long WriteWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,
              char *EndTimeStr,char *StationList);


/*==========================================================================*/
/*              Routines for reading WDC-format files                       */
/*==========================================================================*/


/*--------------------------------------------------------------------------*/
/*  Find size of a data file                                                */
/*--------------------------------------------------------------------------*/

static long GetFileSize(char *FileName)
{
    long    dummy = 0;
    FILE    *DataFile;

    if ((DataFile = fopen(FileName, "r")) != NULL) {
        dummy = fseek(DataFile,0L,SEEK_END);
        dummy = ftell(DataFile);
        fclose(DataFile);
    }
    return dummy;
}


/*--------------------------------------------------------------------------*/
/*  Get the start time in seconds of a 121-character data block. In the     */
/*  WDC file the date is in format YYMMXDD. Where X is the component        */
/*  character (is there any logic here ?).                                  */
/* The two digit century is displayed separately starting at column 14.     */
/*--------------------------------------------------------------------------*/

static Time_sec GetWDCTime(char *Buffer)
{
    char dummyStr[13] = "????????0000";

    // In Nurmijarvi hour WDC files there is no Century but the column 15
    // is used to display the day character, either " ", "S" or "Q".

    if ((Buffer[15] == 'S') || (Buffer[15] == 'Q')) {
        dummyStr[0] = ' ';
        dummyStr[1] = ' ';
    } else {
        strncpy(dummyStr,Buffer+WDC_Century,2);
    }
    strncpy(dummyStr+2,Buffer+WDC_YYMM,4);
    strncpy(dummyStr+6,Buffer+WDC_Day,2);
    return (StrToSecsC(dummyStr,8));
}


/*--------------------------------------------------------------------------*/
/*  Fill some fields in the Station_struct by copying them from the WDC     */
/*  buffer. AddNewStation procedure in the MagnData.h file will set other   */
/*  station parameters.                                                     */
/*--------------------------------------------------------------------------*/

static void SetStationParamsWDC(StationPtr Station,char *Buff)
{
    long dummy;

    /* Set the station ID */
    strncpy(Station->StationID,Buff+WDC_Station,3);
    Station->StationID[3] = '\0';

    /* Set the location of the station */

    if (strncmp("HER",Buff+WDC_Station,3) == 0) {
        Station->Longitude =  1923;
        Station->Latitude  = -3443;
    }
    else
    if (strncmp("BOU",Buff+WDC_Station,3) == 0) {
        Station->Longitude = 25480;
        Station->Latitude  =  4010;
    }
    else
    if (strncmp("HON",Buff+WDC_Station,3) == 0) {
        Station->Longitude = 20200;
        Station->Latitude  =  2130;
    }
    else
    if (strncmp("KAK",Buff+WDC_Station,3) == 0) {
        Station->Longitude = 14018;
        Station->Latitude  =  3623;
    }
    else
    if (strncmp("SJG",Buff+WDC_Station,3) == 0) {
        Station->Longitude = 29380;
        Station->Latitude  =  1810;
    }
    else {
//        fprintf(stderr,"Station longitude not defined\n");
        Station->Longitude = 0;
        Station->Latitude  = 0;
    }
}


/*--------------------------------------------------------------------------*/
/* Get the 4 character number from the Buffer. We could also use sscanf,    */
/* but this is considerably faster because no checkings are made.           */
/*--------------------------------------------------------------------------*/

static long GetWDCFieldValue(char *p)
{
    int k;
    long value = 0;
    long negative = 0;

    for (k=0;k<4;k++) {
        if (*p == '-') negative = 1;
        else
        if (*p != ' ') value = 10*value+(*p-'0');
        p++;
    }

    if (negative) value = -value;
    if (value == MissingWDC) return(MissingValue);
    return (value);
}



/*--------------------------------------------------------------------------*/
/* Get data from the 120 character block into the proper place              */
/*--------------------------------------------------------------------------*/

static void GetWDCData(StationPtr Station, Time_sec Time, char *Buff)
{
    long *DataPtr,*TPtr;
    long i,Count,value,Base;
    char *q;
    char Component;
    char TimeStr[20];

    Component = *(Buff+WDC_Component);

    /* The file may contain H and D as well as X and Y */
    /* so mark the component read                      */
    switch (Component) {
        case 'H' : SetDataValue(Station,Time,'T',1); break;
        case 'X' : SetDataValue(Station,Time,'T',0); break;
        case 'D' : SetDataValue(Station,Time,'A',1); break;
        case 'Y' : SetDataValue(Station,Time,'A',0); break;
    }

    DataPtr = GetDataPtr(Station,Time,Component);

    /* Get the baseline value */
    Base = GetWDCFieldValue(Buff+WDC_Base);
    switch (Component) {
        case 'D' : Base *=  600; break; /* convert into 0.1 arcmin  */
        default  : Base *= 1000; break; /* convert into 0.1 nT      */
    }

    /* Copy the data */
    q = Buff+WDC_Data;
    for (i=0;i<24;i++) {

        value = GetWDCFieldValue(q);

        if (value == MissingValue)
            *DataPtr++ = MissingValue;
        else {
            if (Component == 'D')       /* Unit = 0.1 arcmin */
                *DataPtr++ = Base+value;
            else
                *DataPtr++ = Base+10*value; /* Convert into 0.1 nT */
        }
        q += 4;
    }


    /* Set the component characters in the Station->Components */
    switch (Component) {
        case 'H' :  Station->Components[0] = 'H'; break;
        case 'X' :  Station->Components[0] = 'X'; break;
        case 'D' :  Station->Components[1] = 'D'; break;
        case 'Y' :  Station->Components[1] = 'Y'; break;
        case 'Z' :  Station->Components[2] = 'Z'; break;
    }
}


/*--------------------------------------------------------------------------*/
/* Compute HDF from XYZ and fill the appropriate data fields in the station */
/* structure.                                                               */
/*--------------------------------------------------------------------------*/

void Compute_HDF_from_XYZ(StationPtr s)
{
    long *XPtr,*YPtr,*ZPtr,*HPtr,*DPtr,*FPtr;
    long i,count;
    double X,Y,Z;
    double coeff = 34377.468;       /*  = (180*60*10)/Pi. Unit of D is 0.1 Arcmin */
/*     double coeff = 137509.87; = 4 * (180*60*10)/Pi . Unit of D is (0.1 Arcmin)/4 */

    /* Set the pointers */
    XPtr = s->X;
    YPtr = s->Y;
    ZPtr = s->Z;
    HPtr = s->H;
    DPtr = s->D;
    FPtr = s->F;

    count = (s->EndTime - s->StartTime)/s->TimeStep;
    for (i=0;i<count;i++) {
        X = (double) *XPtr;
        Y = (double) *YPtr;
        Z = (double) *ZPtr;

        if ((*XPtr != MissingValue) && (*YPtr != MissingValue)) {
            *HPtr = RoundFloat(sqrt(X*X+Y*Y));
            *DPtr = RoundFloat(coeff*atan2(Y,X));
        } else {
            *HPtr = MissingValue;
            *DPtr = MissingValue;
        }

        if ((*XPtr != MissingValue) && (*YPtr != MissingValue) && (*ZPtr != MissingValue)) {
            *FPtr = RoundFloat(sqrt(X*X + Y*Y + Z*Z));
        } else {
            *FPtr = MissingValue;
        }

        XPtr++;
        YPtr++;
        ZPtr++;
        HPtr++;
        DPtr++;
        FPtr++;
    }
}

/*--------------------------------------------------------------------------*/
/* Compute XYZ from HDZ and fill the appropriate data fields in the station */
/* structure.                                                               */
/*--------------------------------------------------------------------------*/

void Compute_XYZ_from_HDZ(StationPtr s)
{
    long *XPtr,*YPtr,*ZPtr,*HPtr,*DPtr,*FPtr;
    long i,count;
    double H,D,Z;
    double coeff = 2.90888e-5; /* = Pi/(180*60*10) */
/*  double coeff = 34377.468;         = (180*60*10)/Pi. Unit of D is 0.1 Arcmin */
/*  double coeff = 137509.87; = 4 * (180*60*10)/Pi . Unit of D is (0.1 Arcmin)/4 */

    /* Set the pointers */
    XPtr = s->X;
    YPtr = s->Y;
    ZPtr = s->Z;
    HPtr = s->H;
    DPtr = s->D;
    FPtr = s->F;

    count = (s->EndTime - s->StartTime)/s->TimeStep;
    for (i=0;i<count;i++) {
        H = (double) *HPtr;
        D = (double) *DPtr;
        Z = (double) *ZPtr;

        if ((*HPtr != MissingValue) && (*DPtr != MissingValue)) {
            *XPtr = RoundFloat(H*cos(coeff*D));
            *YPtr = RoundFloat(H*sin(coeff*D));
        } else {
            *XPtr = MissingValue;
            *YPtr = MissingValue;
        }

        if ((*HPtr != MissingValue) && (*ZPtr != MissingValue)) {
            *FPtr = RoundFloat(sqrt(H*H + Z*Z));
        } else {
            *FPtr = MissingValue;
        }

        XPtr++;
        YPtr++;
        ZPtr++;
        HPtr++;
        DPtr++;
        FPtr++;
    }
}


/*--------------------------------------------------------------------------*/
/*  Here is the routine for reading WDC-format data file into memory.       */
/*  The data is read into a Network structure (see MagnData.h for details). */
/*  See the comments at the beginning of this file for more details about   */
/*  the parameters for this routine.                                        */
/*--------------------------------------------------------------------------*/

long ReadWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,
              char *EndTimeStr,char *StationList, long *Century)
{
    FILE    *WDCfile;               /* The data file to be processed        */
    Time_sec CurrTime;              /* Time of the current data block       */
    Time_sec StartTime;             /* Time of first record to be read      */
    Time_sec EndTime;               /* Time of record not read anymore      */
    Time_sec FirstTime,LastTime;    /* Time of record not read anymore      */
    Time_sec T;                     /* Time of the previous block           */
    StationPtr Station;             /* Pointer to the current station       */
    long    WDCFileSize;            /* Size of the WDC file                 */
    char *WDCBuffer;                /* Pointer to the start of the WDC file */
    char *Line;                     /* Pointer to the start of current line */
    char *q;                        /* Dummy pointer                        */
    long DataCount;                 /* Number of data points                */
    long DataSize;                  /* Size of the area allocated to data   */
    char Comp;

    InitNetwork(Network);           /* Initialize the fields of Network     */


    /* ------------------------------------ */
    /*  Try to read the file into memory    */
    /* ------------------------------------ */

    if ((FileName == NULL) || (*FileName == '\0')) {
        fprintf(stderr,"### FileName not specified\n");
        return(FileError);
    }

    WDCFileSize = GetFileSize(FileName);
    if (WDCFileSize == 0) {
        fprintf(stderr,"### Could not open file %s\n",FileName);
        return FileError;
    }

    WDCBuffer = malloc(WDCFileSize+1);
    if (WDCBuffer == NULL) {
        fprintf(stderr,"### Out of memory !\n");
        return(OutOfMemory);
    }
    WDCBuffer[WDCFileSize] = 0; /* Set our end of file mark */

    WDCfile = fopen(FileName,"r");
    if (fread(WDCBuffer,WDCFileSize,1,WDCfile) == 0) {
        fprintf(stderr,"### Could not read file %s\n",FileName);
        return FileError;
    }
    fclose(WDCfile);

    *Century = 10*(WDCBuffer[WDC_Century] - 48) + (WDCBuffer[WDC_Century+1] - 48);


    /* ------------------------------------ */
    /*  Get the times of the first and      */
    /*  last data blocks in the data file.  */
    /* ------------------------------------ */

    FirstTime = GetWDCTime(WDCBuffer);
    q = WDCBuffer+WDCFileSize;  /* q points to the 0 marker at the end  */
    q = q - WDCLineLength;
    if (*q < ' ') q++;      /* check if there are both NL and CR    */
    LastTime  = GetWDCTime(q)+24*3600;      // Time not included anymore


    /* ----------------------------- */
    /*  Set the start and end times  */
    /* ----------------------------- */

    if ((StartTimeStr == NULL) || (*StartTimeStr == '\0'))
        StartTime = FirstTime;
    else StartTime = StrToSecsC(StartTimeStr,0);

    if ((EndTimeStr == NULL) || (*EndTimeStr == '\0'))
        EndTime = LastTime;
    else EndTime = StrToSecsC(EndTimeStr,0);


    /* ---------------------------------- */
    /*  Create the station structure and  */
    /*  set some parameters               */
    /* ---------------------------------- */

    Station = AddNewStation(Network,WDCstep,24*WDCstep);
    SetStationParamsWDC(Station,WDCBuffer);
    FreeOneDayBlock(Station->FirstDay);
    Station->FirstDay  = NULL;
    Station->StartTime = StartTime;
    Station->EndTime   = EndTime;


    /* ---------------------------------- */
    /*  Allocate space for the data and   */
    /*  fill the space with missing data. */
    /* ---------------------------------- */

    DataCount  = (EndTime-StartTime)/WDCstep;
    DataSize   = sizeof(long)*DataCount;
    Station->X = (long *) malloc(DataSize);
    Station->Y = (long *) malloc(DataSize);
    Station->Z = (long *) malloc(DataSize);
    Station->D = (long *) malloc(DataSize);
    Station->H = (long *) malloc(DataSize);
    Station->F = (long *) malloc(DataSize);
    Station->T = (long *) malloc(DataSize/24);
    Station->A = (long *) malloc(DataSize/24);

    if ((Station->X == NULL) || (Station->Y == NULL) ||
        (Station->Z == NULL) || (Station->T == NULL) || (Station->A == NULL))
        return (OutOfMemory);

    FillWithMissingData(Station->X,DataCount);
    FillWithMissingData(Station->Y,DataCount);
    FillWithMissingData(Station->Z,DataCount);
    FillWithMissingData(Station->D,DataCount);
    FillWithMissingData(Station->H,DataCount);
    FillWithMissingData(Station->F,DataCount);
    FillWithMissingData(Station->T,DataCount/24);
    FillWithMissingData(Station->A,DataCount/24);


    /* ------------------------------------------------ */
    /*  Read the data line by line and insert the data  */
    /*  into proper place in the station structure.     */
    /* ------------------------------------------------ */

    Line = WDCBuffer;

    while (*Line != 0)  /* Go through all the lines */
    {
        Comp = Line[WDC_Component];

        if (Comp == '*') Line[WDC_Component] = 'X';  /* for Dst */
        if ((Comp != 'X') && (Comp != 'Y') && (Comp != 'Z') && (Comp != 'H') && (Comp != 'D')  && (Comp != 'F')) return IllegalComponent;

        CurrTime = GetWDCTime(Line);

        if ((CurrTime >= StartTime) && (CurrTime < EndTime) && (Comp != 'F') && (Comp != 'I')) /* skip F and I components */
        {
            GetWDCData(Station,CurrTime,Line);
        }

        Line += WDCLineLength-1;
        if ((*Line > 0) && (*Line < ' ')) Line++; /* If there are both NL and CR */
    }


    free(WDCBuffer);

    if (strcmp(Station->Components,"XYZ") == 0) Compute_HDF_from_XYZ(Station);
    if (strcmp(Station->Components,"HDZ") == 0) Compute_XYZ_from_HDZ(Station);

    return OK;
}


/*==========================================================================*/
/*              Routines for writing WDC-format files                       */
/*==========================================================================*/


/*--------------------------------------------------------------------------*/
/*  Write one 122 character block for given station and for given time      */
/*  into specified datafile.                                                */
/*--------------------------------------------------------------------------*/

static void WriteWDCBlock(FILE *DataFile,StationPtr Station,Time_sec Time,
                          char Component)
{
    long j,Ave;
    char DummyStr[20];
    Time_sec T = Time;
    long MaxValue,MinValue;
    long BaseValue;

    /*** Write station info and time fields  ***/
    fprintf(DataFile,"%s",Station->StationID);
    SecsToStrC(Time,DummyStr);  // YYYYMMDD
    fprintf(DataFile,"%.4s",DummyStr+2);
    fprintf(DataFile,"%c",Component);
    fprintf(DataFile,"%.2s",DummyStr+4);
    fprintf(DataFile,"      "); /* 6 spaces */

    /*** Get the baseline value ***/
    FindMaxMin(Station,Component,T,24*3600, &MaxValue, &MinValue);

    if (MinValue == MissingValue)
        BaseValue = MissingValue;
    else {
        if (Component == 'D')
            BaseValue = MinValue/600;   /* Convert into degrees  */
        else
            BaseValue = MinValue/1000;  /* Convert into 100 nT's */
    }

    if (BaseValue == MissingValue)
        fprintf(DataFile,"9999");
    else
        fprintf(DataFile,"% 4d",BaseValue);

    /*** Write field values ***/
    for (j=0;j<24;j++) {
        Ave = GetDataValue(Station,T,Component);
        if (Ave == MissingValue)
            fprintf(DataFile,"9999");
        else {
            if (Component == 'D')       /* Unit = 0.1 arcmin */
                fprintf(DataFile,"% 4d",Ave-600*BaseValue);
            else
                fprintf(DataFile,"% 4d",Ave/10 - 100*BaseValue);    /* Unit = 1 nT */
        }
        T += 3600;
    }

    /*** Write hour average ***/
    Ave = ComputeAverage(Station,Component,Time,24*3600);
    if (Ave != MissingValue) {
        if (Component == 'D')       /* Unit = 0.1 arcmin */
            fprintf(DataFile,"% 4d",Ave-600*BaseValue);
        else
            fprintf(DataFile,"% 4d",Ave/10 - 100*BaseValue);
    }
    else
        fprintf(DataFile,"9999");

     /*** The block end characters ***/
    fprintf(DataFile,"\n");
}


/*--------------------------------------------------------------------------*/
/*  Here is the routine for writing data into an WDC-format file. See the   */
/*  comments at the beginning of this file for more details about the       */
/*  parameters for this routine.                                            */
/*--------------------------------------------------------------------------*/

long WriteWDChour(char *FileName,NetworkPtr Network,char *StartTimeStr,
              char *EndTimeStr,char *StationList)
{
    FILE    *WDCfile;               /* The data file to be processed        */
    Time_sec Time;                  /* Time of the current data block       */
    Time_sec StartTime;             /* Time of first record to be read      */
    Time_sec EndTime;               /* Time of record not read anymore      */
    long    Comp;                   /* Dummy index                          */
    StationPtr Station;             /* Pointer to the current station       */
    Time_struct TimeStruct;

    /* --- Check there is some data in the Network --- */
    if (Network->StationList == NULL) return FAIL;

    /* --- Set the start and end times --- */
    if ((StartTimeStr == NULL) || (*StartTimeStr == '\0'))
        StartTime = Network->StationList[0].StartTime;
    else StartTime = StrToSecsC(StartTimeStr,0);

    if ((EndTimeStr == NULL) || (*EndTimeStr == '\0'))
        EndTime = Network->StationList[0].EndTime;
    else EndTime = StrToSecsC(EndTimeStr,0);

    /* --- Try to create the file --- */
    if ((FileName == NULL) || (*FileName == '\0')) WDCfile = stdout;
    else if ((WDCfile = fopen(FileName, "w")) == NULL) return FileError;

    /* ------------------- Write the data ------------------*/
    /* Data is written one day per line (24 hour values).   */
    /* -----------------------------------------------------*/

    Time = StartTime;
    while (Time < EndTime) {
        Station = Network->StationList;
        while (Station != NULL) {
            if (StationInList(Station->StationID,StationList)) {
                for (Comp=0;Comp<3;Comp++) {
                    WriteWDCBlock(WDCfile,Station,Time,Station->Components[Comp]);
                }
            }
            Station = Station->Next;
        }
        Time += 24*3600;
    }
    fclose(WDCfile);

    return OK;
}