/****************************************************************************/
/*                                                                          */
/*                                  NOR.h                                   */
/*                                                                          */
/****************************************************************************/
/****************************************************************************/
/* This is a C language header file that defines the following routines for */
/* reading data from a NOR-format data file into memory.					*/
/* The NOR format the new data format used by Auroral Observatory of the	*/
/* University of NORmsa since January 2009 in storing magnetometer data		*/
/* from their stations (e.g. Svalbard stations + danish stations). The name */
/* NOR is not the official name of this data format but since I haven't     */
/* seen any name used for this data format I will call it NOR (=Norway).    */
/* These routines will handle correctly 1 second (16 bit), 10 second		*/
/* and one minute data files but not 1 second (32 bit) files.				*/
/*                                                                          */
/*  -------------------  Read data from a NOR-format file   --------------  */
/*                                                                          */
/* long ReadNOR(char *FileName, NetworkPtr Network, char *StartTimeStr,     */
/*               char *EndTimeStr, char *StationList, Time_sec SampleStep)  */
/*      FileName        Name of the NOR-file containing the data.           */
/*                      If FileName is nil or zero length then standard     */
/*                      input is used.                                      */
/*      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 included.              */
/*		SampleStep		Time separation between successive data points.		*/
/*                                                                          */
/*      The function result indicates how successfull the file reading was: */
/*          0: OK               : File read successfully                    */
/*          1: FileError        : Failed to read the given file             */
/*          2: OutOfMemory      : Memory allocation failed during reading   */
/*          3: FileFormatError  : Data file is not a NOR file               */
/*                                                                          */
/*                                                                          */
/*  ----------------------------------------------------------------------- */
/*                                                                          */
/*  NOTE:   These routines are intended to be used only in dealing with     */
/*          IMAGE magnetometer network data files in NOR-format.            */
/*          The routines defined here make assumptions about the data       */
/*          content. The routines might work also with other than IMAGE     */
/*          data files but this is not guaranteed.                          */
/*                                                                          */
/****************************************************************************/
/****************************************************************************/
/*                            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.0  29.7.2009                              */
/****************************************************************************/
/****************************************************************************/
/*  Version history:                                                        */
/*                                                                          */
/*  1.0  29.7.2009  First official release                                  */
/****************************************************************************/


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


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

// Offsets of the various subblocks
#define NOR_Hd_block		 0		// Offset of the header block
#define NOR_S1_block	   117		// Offset of the scale values #1 block
#define NOR_Q1_block	   249		// Offset of the QMN values #1 block
#define NOR_S2_block	   315		// Offset of the scale values #2 block
#define NOR_Q2_block	   447		// Offset of the QMN values #2 block
#define NOR_data_block	   513		// Offset of the data block

// Offsets of the various fields within header block
#define NOR_StationNbr       0
#define NOR_StationName      1
#define NOR_StationID	    22
#define NOR_Latitude     	28
#define NOR_Longitude     	39
#define NOR_Height          50
#define NOR_Year            52
#define NOR_Month           54
#define NOR_Day             56
#define NOR_Hour            58
#define NOR_Temp_magn       60
#define NOR_Temp_sensor     66
#define NOR_Exitation       72
#define NOR_Origin          78
#define NOR_Timing_Origin   79
#define NOR_Instrument      80
#define NOR_Filter         101
#define NOR_Timestamp      102
#define NOR_Free	       103
#define NOR_Dno1           113
#define NOR_Dno2           115

// Offsets of the various fields within Scale values block
#define NOR_S_Xa             0
#define NOR_S_Ya            11
#define NOR_S_Za            22
#define NOR_S_X0            33
#define NOR_S_Y0            44
#define NOR_S_Z0            55
#define NOR_S_D0            66
#define NOR_S_I0            77
#define NOR_S_N1            88
#define NOR_S_N2            99
#define NOR_S_N3           110
#define NOR_S_N4           121

// Offsets of the various fields within QMNvalues block
#define NOR_Q_QmX            0
#define NOR_Q_QmY           11
#define NOR_Q_QmZ           22
#define NOR_Q_QmD           33
#define NOR_Q_QmI           44
#define NOR_Q_QmF           55

#define MissingNOR		-32768
#define Deg_to_Rad		0.0174532925

typedef char *NORptr;


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

static void P_to_C_str(char *P_str, char *C_str);
static double atof_P_str(char *P_str);
static long GetWordLH(char *p);
static long GetShortLH(char *p);
static Time_sec GetNORTime(NORptr NOR);
static long ReadNORBlock(FILE *NORfile, NORptr NOR, long BlockSize);
static void GetNORStationID(char *StationID,NORptr NOR);
static void SetStationParamsNOR(StationPtr Station,NORptr NOR, Time_sec StartTime, long SampleStep);
static void FindScaleValues(NORptr Block, double *Xscale, double *Yscale, double *Zscale,
							double *D0, double *I0, double *cosI0, double *sinI0);
static void GetNORData(StationPtr Station,NORptr NOR);
long ReadNOR(char *FileName, NetworkPtr Network, char *StartTimeStr,
              char *EndTimeStr, char *StationList, Time_sec SampleStep);


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

/*--------------------------------------------------------------------------*/
/*  Convert a Pascal string into a C string									*/
/*--------------------------------------------------------------------------*/

static void P_to_C_str(char *P_str, char *C_str)
{
    register long i,count;
    register char *p = P_str+1;
    register char *q = C_str;

    count = *P_str;
    for (i=0;i<count;i++) *q++ = *p++;
    *q = '\0';
}

/*--------------------------------------------------------------------------*/
/*  Convert a Pascal string representing a floating point number into double*/
/*--------------------------------------------------------------------------*/

static double atof_P_str(char *P_str)
{
    char C_str[20];

    P_to_C_str(P_str,C_str);
    return (atof(C_str));
}


/*--------------------------------------------------------------------------*/
/*  Get a two-byte binary number. The order of the bytes is : Low-High      */
/*  We cannot read the binaries with (short) p, where p is defined as       */
/*  (short) *p, since some processors use high-low byte ordering. In the    */
/*  NOR format the ordering is defined as low-high. We could check the      */
/*  type of the processor and then swap the bytes if necessary but that     */
/*  wouldn't be any simpler.                                                */
/*--------------------------------------------------------------------------*/

static long GetWordLH(char *p)		// Get unsigned integer
{
    register long high,low;

    low  = 0x000000FF & (*p++);
    high = 0x000000FF & (*p);
    return((high << 8) | low);
}

static long GetShortLH(char *p)		// Get signed integer
{
    register long high,low;

    low  = 0x000000FF & (*p++);
    high = (signed char) (*p);
    return((high << 8) | low);
}


/*--------------------------------------------------------------------------*/
/*  Get the start time in seconds of a NOR data block.                      */
/*--------------------------------------------------------------------------*/

static Time_sec GetNORTime(NORptr NOR)
{
    Time_struct MyTime;

    MyTime.tm_year  = GetWordLH(NOR+NOR_Year)-1900;
    MyTime.tm_mon   = GetWordLH(NOR+NOR_Month)-1;
    MyTime.tm_mday  = GetWordLH(NOR+NOR_Day);
    MyTime.tm_hour  = GetWordLH(NOR+NOR_Hour);
    MyTime.tm_min   = 0;
    MyTime.tm_sec   = 0;

// fprintf(stderr,"%04d %02d %02d   %02d\n",MyTime.tm_year+1900,MyTime.tm_mon+1,MyTime.tm_mday,MyTime.tm_hour);

    return TmToSecs(&MyTime);
}


/*--------------------------------------------------------------------------*/
/* Read one data block from the NOR-format file into NOR block. Function 	*/
/* returns 1 if reading was successful and 0 if the end of file is reached. */
/*--------------------------------------------------------------------------*/

static long ReadNORBlock(FILE *NORfile, NORptr NOR, long BlockSize)
{
    long FileStatus;

    FileStatus = fread(NOR,BlockSize,1,NORfile);
    if (FileStatus == 0)
    	return (0);
    else
    	return (1);
}


/*--------------------------------------------------------------------------*/
/*  Find out the three letter station ID from the ID-number code in the     */
/*  NOR block.                                                              */
/*--------------------------------------------------------------------------*/

static void GetNORStationID(char *StationID,NORptr NOR)
{
    char ID[10];
    int i;

//	Before 2011 the station ID was determined from one byte station number.
//	With inclusion of Danish magnetometer data the system changed and the
//	station number was put to 99 (= missing station) for all stations.
//     ID = *(NOR+NOR_Hd_block+NOR_StationNbr);
//     while (i<IMAGEStationCount) {
//         if (StationInfoTable[i++].TRO_ID == ID) break;
//     }
//     strcpy(StationID,StationInfoTable[i-1].StationID);

	P_to_C_str(NOR+NOR_Hd_block+NOR_StationID,ID);
	for (i=0;i<3;i++) StationID[i] = toupper(ID[i]);
	StationID[3] = '\0';
}


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

static void SetStationParamsNOR(StationPtr Station,NORptr NOR, Time_sec StartTime, long SampleStep)
{
    char Lat_str[10],Lon_str[10];

    /* --- Set the three letter station ID --- */
    GetNORStationID(Station->StationID,NOR);
    Station->StationID[3] = '\0';

    /* --- Set the coordinates of the station, unit = 0.01 degrees --- */
    P_to_C_str(NOR+NOR_Latitude,Lat_str);
    P_to_C_str(NOR+NOR_Longitude,Lon_str);

    Station->Latitude = RoundFloat(100*atof(Lat_str));
    Station->Longitude = RoundFloat(100*atof(Lon_str));

	/* Set the component markers */
	strcpy(Station->Components,"XYZ");

    /* --- Set StartTime and TimeStep, GetNORData-routine sets EndTime --- */
    Station->StartTime = StartTime;
    Station->EndTime   = Station->StartTime;
    Station->TimeStep  = SampleStep;
}


/*--------------------------------------------------------------------------*/
/*  Find the scale values , offsets and declination and inclination values  */
/*	from the Scale values data block.										*/
/*--------------------------------------------------------------------------*/

static void FindScaleValues(NORptr Block, double *Xscale, double *Yscale, double *Zscale,
							double *D0, double *I0, double *cosI0, double *sinI0)
{
    Xscale[0] = atof_P_str(Block + NOR_S_Xa);	// index 0: scale values
    Yscale[0] = atof_P_str(Block + NOR_S_Ya);
    Zscale[0] = atof_P_str(Block + NOR_S_Za);
    Xscale[1] = atof_P_str(Block + NOR_S_X0);	// index 1: offsets
    Yscale[1] = atof_P_str(Block + NOR_S_Y0);
    Zscale[1] = atof_P_str(Block + NOR_S_Z0);
	*D0 = atof_P_str(Block + NOR_S_D0);
	*I0 = atof_P_str(Block + NOR_S_I0);
	*sinI0 = sin((*I0)*Deg_to_Rad);				// D0 and I0 are in degrees
	*cosI0 = cos((*I0)*Deg_to_Rad);
}


/*--------------------------------------------------------------------------*/
/*  Copy data from NOR block into proper OneDayBlock in Station structure.  */
/*--------------------------------------------------------------------------*/

static void GetNORData(StationPtr Station,NORptr NOR)
{
    OneDayPtr p;
    long *Xptr, *Yptr, *Zptr, *Tptr;
	double Xscale[2];	// scales and offsets. 0 = scale, 1 = offset
	double Yscale[2];
	double Zscale[2];
	double D0, I0, cosI0, sinI0;
	double XX, YY, ZZ, tx, H, D;
	long Scale_2_start;
	long Count;
	long DataOrigin;
	long i, step;
	long X,Y,Z;

	/* Find the last OneDayBlock */
	for (p = Station->FirstDay; p->NextDay != NULL; p = p->NextDay);

    /* Set the temperature value */
	Tptr  = (p->T)+(p->TCount);
	*Tptr = GetShortLH(NOR + NOR_Temp_magn);
	if (*Tptr == MissingNOR) *Tptr = MissingValue;
	p->TCount++;

    /* Get the baselines and scalefactors for the first block */
	FindScaleValues(NOR+NOR_S1_block, Xscale, Yscale, Zscale, &D0, &I0, &cosI0, &sinI0);

	Scale_2_start = GetWordLH(NOR + NOR_Dno1);	// Start index of second scale value set


    /* Get the number of datapoints */
    switch (Station->TimeStep) {
    	case 1  : Count = 10800; break;	// 1 second short
    	case 10 : Count =  1080; break;	// 10 seconds
    	case 60 : Count =  1440; break;	// 1 minute
		default : Count =  1080; 		// 10 seconds
	}

    /* Get raw data, multiply by scale coefficients, add offsets */
    /* and rotate into geographic coordinates if necessary.		 */

    DataOrigin = *(NOR + NOR_Origin);	// 0 = raw data, 1 = geographic XYZ

    Xptr = (p->X)+(p->XCount);
    Yptr = (p->Y)+(p->YCount);
    Zptr = (p->Z)+(p->ZCount);
    step = 0;

    for (i=0;i<Count;i++) {
        X = GetShortLH(NOR + NOR_data_block + step);
        Y = GetShortLH(NOR + NOR_data_block + 2*Count + step);
        Z = GetShortLH(NOR + NOR_data_block + 4*Count + step);

        if ((X == MissingNOR) || (Y == MissingNOR) || (Z == MissingNOR)) {
        	*Xptr++ = MissingValue;
        	*Yptr++ = MissingValue;
        	*Zptr++ = MissingValue;
        }
        else {
        	XX = Xscale[0]*X + Xscale[1];
        	YY = Yscale[0]*Y + Yscale[1];
        	ZZ = Zscale[0]*Z + Zscale[1];

			if (DataOrigin == 0) {		// Raw data, rotate to geographic coordinates
				tx = ZZ*cosI0 + XX*sinI0;			// temporary 'X'
				H  = sqrt(tx*tx + YY*YY);
				D  = atan2(YY,tx) + Deg_to_Rad*D0;	// Declination in radians
				*Xptr++ = RoundFloat(10.0*H*cos(D));				// X geographic
				*Yptr++ = RoundFloat(10.0*H*sin(D));				// Y geographic
				*Zptr++ = RoundFloat(10.0*(ZZ*sinI0 - XX*cosI0));	// Z geographic
			}
			else {						// already geographic XYZ
				*Xptr++ = RoundFloat(10.0*XX);		// Unit = 0.1 nT
				*Yptr++ = RoundFloat(10.0*YY);
				*Zptr++ = RoundFloat(10.0*ZZ);
			}
        }

		if (i == (Scale_2_start-1))		// Check if we need to update the set of coefficients
			FindScaleValues(NOR+NOR_S2_block, Xscale, Yscale, Zscale, &D0, &I0, &cosI0, &sinI0);

        step += 2;
    }

    /* --- Update counters --- */
	p->XCount += Count;
	p->YCount += Count;
	p->ZCount += Count;

    /* --- Set the end time --- */
    Station->EndTime = GetNORTime(NOR)+Count*Station->TimeStep;
}


/*--------------------------------------------------------------------------*/
/* Add missing data block 													*/
/*--------------------------------------------------------------------------*/

static void AddMissingNORBlock(StationPtr Station, long Count)
{
	OneDay_struct *p;

	/* Find the last OneDayBlock */
	for (p = Station->FirstDay; p->NextDay != NULL; p = p->NextDay);

	/* We only have to update counters since the data areas are already */
	/* initialized with missing data									*/
	p->XCount += Count;
	p->YCount += Count;
	p->ZCount += Count;
	p->TCount++;
}



/*--------------------------------------------------------------------------*/
/*  Here is the routine for reading NOR-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.                                        */
/*                                                                          */
/*  The rounite works by reading the data file one block at a time and adds */
/*  that data in the particular stations data block.					    */
/*  For each new station that is encountered space for one day is allocated.*/
/*  If the day becomes full a new one day block is allocated. At the end of */
/*  the routine all one day blocks are combined into a single data block.   */
/*--------------------------------------------------------------------------*/

long ReadNOR(char *FileName, NetworkPtr Network, char *StartTimeStr,
              char *EndTimeStr, char *StationList, Time_sec SampleStep)
{
    FILE    *NORfile;               /* The data file to be processed        */
    char    *NORblock;         		/* Pointer to data record			    */
    Time_sec CurrTime;              /* Time of the current data block       */
    Time_sec PrevTime = MIN_TIME;   /* Time of the previous data block      */
    Time_sec StartTime;             /* Time of first record to be read      */
    Time_sec EndTime;               /* Time of record not read anymore      */
    StationPtr Station;             /* Pointer to the current station       */
    char    StationID[4];           /* Station ID of the current station    */
    long    Count;                  /* Number of data points in one block   */
    long    BlockStatus;            /* Status number for block reading      */
    long	BlockSize;				/* Size of a single data block			*/
char CurrTimeStr[20];
    InitNetwork(Network);           /* Initialize the fields of Network     */

    /* --- Set the start and end times --- */
    if ((StartTimeStr == NULL) || (StartTimeStr[0] == '\0'))
        StartTime = MIN_TIME;
    else StartTime = StrToSecs(StartTimeStr,0);

    if ((EndTimeStr == NULL) || (EndTimeStr[0] == '\0'))
        EndTime = MAX_TIME-(Time_sec)86400;
    else EndTime = StrToSecs(EndTimeStr,0);

    /* --- Allocate space for data record --- */
    switch (SampleStep) {
    	case 1  : BlockSize = NOR_data_block + 64800; Count = 10800; break;	// 1 second short
    	case 10 : BlockSize = NOR_data_block + 6480;  Count =  1080; break;	// 10 seconds
    	case 60 : BlockSize = NOR_data_block + 8640;  Count =  1440; break;	// 1 minute
		default : BlockSize = NOR_data_block + 6480;  Count =  1080; 		// 10 seconds
	}

	NORblock = (char *) malloc(BlockSize);


    /* --- Try to open the file --- */
    if ((FileName == NULL) || (FileName[0] == '\0')) NORfile = stdin;
    else if ((NORfile = fopen(FileName, "rb")) == NULL) return(FileError);

    /* --- Read data into one day blocks ---*/
    while (BlockStatus = ReadNORBlock(NORfile,NORblock,BlockSize))
    {
        if ((CurrTime=GetNORTime(NORblock)) >= EndTime+(Time_sec)86400) break;
// SecsToStr(CurrTime,CurrTimeStr);
// fprintf(stderr,"%s\n",CurrTimeStr);


        GetNORStationID(StationID,NORblock);

        if ((CurrTime >= PrevTime) && (CurrTime >= StartTime) && (CurrTime < EndTime)
            && StationInList(StationID,StationList))
        {
            Station = FindStation(Network,StationID);
            if (Station == NULL) {      /* First occurrence of this station */
                 Station = AddNewStation(Network,SampleStep,Count*SampleStep);
                if (Station != NULL) {
                    SetStationParamsNOR(Station,NORblock,CurrTime,SampleStep);
				}
                else        /* Unable to allocate memory for station data   */
                    return(OutOfMemory);
            }

			/* --- Add possible missing data blocks ---                */
			/* PrevTime is the expected time of the current data block.*/
			/* If it is less than the time read from the block then    */
			/* add missing data blocks                                 */

			PrevTime = Station->EndTime;
			while (PrevTime < CurrTime) {
				AddMissingNORBlock(Station, Count);
				if (OneDayFull(Station,'Z')) {
					if (AddOneDay(Station))
						return(OutOfMemory);
				}
				PrevTime += Count*SampleStep;
			}

            if (OneDayFull(Station,'X')) {
                if (AddOneDay(Station))
                    return(OutOfMemory);    /* Unable to allocate memory    */
            }
            GetNORData(Station,NORblock);
        }
    }
    fclose(NORfile);

    /* --- Combine one day blocks into one large block ---*/
    Station = Network->StationList;
    while (Station != NULL) {
        if (CombineData(Station) == OutOfMemory) return(OutOfMemory);
        Station = Station->Next;
    }

    return OK;
}