/****************************************************************************/
/*                                                                          */
/*                           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] [<] IAGA_file			*/
/*         [-s YYMMDD]          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 YYMMDD]          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.               */
/*         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.0        15.6.2010                        */
/****************************************************************************/
/****************************************************************************/
/*  Version history:                                                        */
/*                                                                          */
/*  1.0  15.06.2010 First official release                                  */
/****************************************************************************/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include "Usage.h"
#include "MagnData.h"
#include "K_index.h"
#include "IAGA.h"

#define OneDay 86400

char *version = "1.0";
char *date = "15.6.2010";


#define HeaderLineCount 5
#define UsageLineCount  19

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 YYMMDD]       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 YYMMDD]       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                       ",
"     IAGA_file         Name of IAF-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);
    }
}



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

void Compute_HDF(StationPtr s)
{
    long *XPtr,*YPtr,*ZPtr,*HPtr,*DPtr,*FPtr;
    long i,count;
    double X,Y,Z,H;

    /* 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)) {
        	H = sqrt(X*X+Y*Y);
            *HPtr = RoundFloat(H);
            *DPtr = RoundFloat(H*atan2(Y,X));	// D(nT) = H*D(rad)
        } 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++;
    }
}



/*--------------------------------------------------------------------------*/
/*                          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[14] = "";
    char	EndTimeStr[14]   = "";
    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 '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;
    }

    if (strlen(StartTimeStr) == 0) {
		Station = NETWORK.StationList;
		SecsToStr(Station->StartTime + OneDay,StartTimeStr);
		SecsToStr(Station->EndTime - 2*OneDay,EndTimeStr);
    }
	else if (strlen(EndTimeStr) == 0) strcpy(EndTimeStr,StartTimeStr);



    /*=======================================*/
    /* Go through all the days               */
    /*=======================================*/

	StartTime = StrToSecs(StartTimeStr,0);
	EndTime = StrToSecs(EndTimeStr,0);

	Station = NETWORK.StationList;
	if ((StartTime < Station->StartTime + OneDay) || (EndTime > Station->EndTime - 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 */
            SecsToStr(T,TimeStr); TimeStr[6] = '\0';
			if (T >= YEAR2000)
				printf("20");
			else
				printf("19");
            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;
}