/*
 * EasyWave - A realtime tsunami simulation program with GPU support.
 * Copyright (C) 2014  Andrey Babeyko, Johannes Spazier
 * GFZ German Research Centre for Geosciences (http://www.gfz-potsdam.de)
 *
 * Parts of this program (especially the GPU extension) were developed
 * within the context of the following publicly funded project:
 * - TRIDEC, EU 7th Framework Programme, Grant Agreement 258723
 *   (http://www.tridec-online.eu)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <stdio.h>
#include <ctype.h>
#include <string.h>
#include <utilits.h>
#include <cSphere.h>


//=========================================================================
// Constructor
cObsArray::cObsArray()
{
  nPos = 0;
  nObs = 0;
  id = NULL;
  lon = NULL;
  lat = NULL;
  obs = NULL;
}


//=========================================================================
// Destructor
cObsArray::~cObsArray()
{
  if( lon ) delete [] lon;
  if( lat ) delete [] lat;
  if( obs ) {
    for( int n=0; n<nPos; n++ )
      delete [] obs[n];
    delete [] obs;
  }
  if( id ) {
    for( int n=0; n<nPos; n++ )
      delete [] id[n];
    delete [] id;
  }

}


//=========================================================================
// Reset observation array
int cObsArray::resetObs()
{
  if( !obs ) return 0;

  for( int n=0; n<nPos; n++ ) {
    memset( obs[n], 0, nObs*sizeof(double) );
  }

  return 0;
}


//=========================================================================
// Fully reset observation array to a new number of observations
int cObsArray::resetObs( int newnobs )
{

  if( obs ) {
    for( int n=0; n<nPos; n++ )
      delete [] obs[n];
    delete [] obs;
  }

  nObs = newnobs;
  obs = new double* [nPos];
  for( int n=0; n<nPos; n++ ) {
    obs[n] = new double [nObs];
    memset( obs[n], 0, nObs*sizeof(double) );
  }

  return 0;
}


//=========================================================================
// Find poi
int cObsArray::findById( char *id0 )
{
  int n;

  for( n=0; n<nPos; n++ )
    if( !strcmp( id0, id[n] ) ) return n;

  return -1;
}


//=========================================================================
// Read observations from text file
int cObsArray::read( char *fname )
{
  FILE *fp;
  char record[256], buf[64];
  int n,k,idPresent,line=0;
  double lonr,latr;


  fp = fopen(fname,"rt");
  if( fp == NULL ) return Err.post(Err.msgOpenFile(fname));

  // Get number of values per location
  line = 0;
  if( utlReadNextRecord( fp, record,  &line ) == EOF ) return Err.post("Unexpected EOF: %s", fname);
  fclose(fp);
  // check if site ID's are available in the first column. Criterium: first character is not digit
  sscanf( record, "%s", buf );
  if( isalpha( buf[0] ) )
    idPresent = 1;
  else
    idPresent = 0;
  // nObs = number of columns minus 2 (lon lat) minus 1 (if idPresent)
  nObs = utlCountWordsInString( record ) - 2 - idPresent;
  if( nObs < 0 ) return Err.post(Err.msgReadFile( fname, line, "expected: [id lon lat [obs...]" ));

  // Get number of positions
  nPos = utlGetNumberOfRecords( fname );

  // Allocate memory
  lon = new double [nPos];
  lat = new double [nPos];
  if( idPresent ) {
    id = new char* [nPos];
    for( n=0; n<nPos; n++ )
      id[n] = NULL;
  }
  if( nObs > 0 ) {
    obs = new double* [nPos];
    for( n=0; n<nPos; n++ ) {
      obs[n] = new double [nObs];
      memset( obs[n], 0, nObs*sizeof(double) );
    }
  }

  // Read data
  fp = fopen(fname,"rt");
  line = 0;
  for( n=0; n<nPos; n++ ) {
    if( utlReadNextRecord( fp, record,  &line ) == EOF ) return Err.post("Unexpected EOF: %s", fname);

    if( utlCountWordsInString( record ) != (2+idPresent+nObs) ) return Err.post(Err.msgReadFile( fname, line, "invalid number of values" ));

    if( idPresent ) {
      if( sscanf( record, "%s %lf %lf", buf, &lon[n], &lat[n] ) != 3 ) return Err.post(Err.msgReadFile( fname, line, "expected: id lon lat obs..." ));
      id[n] = strdup(buf);
    }
    else {
      if( sscanf( record, "%lf %lf", &lon[n], &lat[n] ) != 2 ) return Err.post(Err.msgReadFile( fname, line, "expected: lon lat obs..." ));
    }

    for( k=0; k<nObs; k++ ) {
      if( utlPickUpWordFromString( record, 3+idPresent+k, buf ) != 0 ) return Err.post(Err.msgReadFile( fname, line, "expected: id lon lat obs..." ));
      if( sscanf( buf, "%lf", &obs[n][k] ) != 1 ) return Err.post(Err.msgReadFile( fname, line, "expected: id lon lat obs..." ));
    }
  }
  fclose( fp );

  return 0;
}


//=========================================================================
// Write to simple text file
int cObsArray::write( char *fname )
{
  FILE *fp;
  int n,k;


  fp = fopen(fname,"wt");

  for( n=0; n<nPos; n++ ) {

    if( id )
      fprintf( fp, "%s %g %g", id[n], lon[n],lat[n] );
    else
      fprintf( fp, "%g %g", lon[n],lat[n] );

    for( k=0; k<nObs; k++ )
      fprintf( fp, " %g", obs[n][k] );

    fprintf( fp, "\n" );
  }

  fclose( fp );

  return 0;
}


//=========================================================================
// Write to binary stream
long cObsArray::writeBin( FILE *fp )
{
  long bytes_written;
  float fbuf;


  bytes_written = 0;
  for( int n=0; n<nPos; n++ ) {
    for( int k=0; k<nObs; k++ ) {
      fbuf = (float)obs[n][k];
      fwrite( &fbuf, sizeof(float), 1, fp );
      bytes_written += sizeof(float);
    }
  }

  return bytes_written;
}



//=========================================================================
// Read from binary stream
long cObsArray::readBin( FILE *fp )
{
  long bytes_read;
  float fbuf;


  bytes_read = 0;
  for( int n=0; n<nPos; n++ ) {
    for( int k=0; k<nObs; k++ ) {
      if( fread( &fbuf, sizeof(float), 1, fp ) != 1 )
        return utlPostError("Unexpected EOF");
      obs[n][k] = (double)fbuf;
      bytes_read += sizeof(float);
    }
  }

  return bytes_read;
}



//=========================================================================
// Calculate observation residual
double cObsArray::residual( cObsArray& ref )
{
  double resid=0.;

  for( int n=0; n<nPos; n++ ) {
    for( int k=0; k<nObs; k++ ) {
      resid += pow( (obs[n][k] - ref.obs[n][k]), 2. );
    }
  }

  resid = sqrt( resid )/nPos/nObs;

  return resid;
}



//=========================================================================
// Calculate norm of observations
double cObsArray::norm()
{
  double norm=0.;

  for( int n=0; n<nPos; n++ ) {
    for( int k=0; k<nObs; k++ ) {
      norm += obs[n][k]*obs[n][k];
    }
  }

  norm = sqrt( norm )/nPos/nObs;

  return norm;
}


// Haversine formula for distance between any two points on the Earth surface
double GeoDistOnSphere( const double lon1, const double lat1, const double lon2, const double lat2 )
{
  const double G2R = 3.14159265358979/180.;  // multiplyer to convert from degrees into radians
  const double REARTH = 6378.137;  // Earth radius in km along equator
  double a,c,dist,rad;

  a = pow( sin(G2R*(lat2-lat1)/2), 2. ) + cos(G2R*lat1)*cos(G2R*lat2)*pow( sin(G2R*(lon2-lon1)/2), 2. );
  rad = sqrt(a);
  if( rad > 1 ) rad = 1;
  c = 2 * asin(rad);
  dist = REARTH*c;

  return dist;
}

double GeoStrikeOnSphere( double lon1, double lat1, double lon2, double lat2 )
{
  const double G2R = 3.14159265358979/180.;  // multiplyer to convert from degrees into radians
  const double R2G = 180./3.14159265358979;  // multiplyer to convert from radians into degrees
  const double REARTH = 6378.137;  // Earth radius in km along equator
  double strike, distRad;


  if( (lat1 == lat2) && (lon1 == lon2) ) {
    strike = 0.;
  }
  else if( lon1 == lon2 ) {
    if( lat1 > lat2 )
      strike = 180.;
    else
      strike = 0.;
  }
  else {
    distRad = GeoDistOnSphere( lon1,lat1,lon2,lat2 ) / REARTH;
    strike = R2G * asin( cos(G2R*lat2)*sin(G2R*(lon2-lon1)) / sin(distRad) );

    if( (lat2 > lat1) && (lon2 > lon1) ) {
    }
    else if( (lat2 < lat1) && (lon2 < lon1) ) {
      strike = 180.0 - strike;
    }
    else if( (lat2 < lat1) && (lon2 > lon1) ) {
      strike = 180.0 - strike;
    }
    else if( (lat2 > lat1) && (lon2 < lon1) ) {
      strike += 360.0;
    }

  }

//  if( strike > 180.0 ) strike -= 180.0;

  return strike;
}