/*
 * 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 <string.h>
#include <utilits.h>
#include <cOgrd.h>
#include "cOkadaEarthquake.h"
#include "easywave.h"

int Imin;
int Imax;
int Jmin;
int Jmax;

#define SRC_GRD 1
#define SRC_FLT 2

//====================================================
int ewSource()
{
  char dsaa_label[8];
  int i,j,ierr,srcType;
  double lon,lat,dz,absuzmax,absuzmin;
  FILE *fp;
  cOkadaEarthquake eq;
  cOgrd uZ;

  CNode& Node = *gNode;

  // check input file type: GRD or fault
  if( (fp = fopen( Par.fileSource, "rb" )) == NULL ) return Err.post( Err.msgOpenFile(Par.fileSource) );
  memset( dsaa_label, 0, 5 );
  ierr = fread( dsaa_label, 4, 1, fp );
  if( !strcmp( dsaa_label,"DSAA" ) || !strcmp( dsaa_label,"DSBB" ) )
    srcType = SRC_GRD;
  else
    srcType = SRC_FLT;
  fclose(fp);


  // load GRD file
  if( srcType == SRC_GRD) {
    ierr = uZ.readGRD( Par.fileSource ); if(ierr) return ierr;
  }

  // read fault(s) from file
  if( srcType == SRC_FLT) {
    int effSymSource = 0;
    long l;
    double dist,energy,factLat,effRad,effMax;
    
    ierr = eq.read( Par.fileSource ); if(ierr) return ierr;

    if( Par.adjustZtop ) {

		// check fault parameters
		Err.disable();
		ierr = eq.finalizeInput();
		while( ierr ) {
		  i = ierr/10;
		  ierr = ierr - 10*i;
		  if( ierr == FLT_ERR_STRIKE ) {
			Log.print( "No strike on input: Employing effective symmetric source model" );
			if( eq.nfault > 1 ) { Err.enable(); return Err.post("Symmetric source assumes only 1 fault"); }
			eq.fault[0].strike = 0.;
			effSymSource = 1;
		  }
		  else if( ierr == FLT_ERR_ZTOP ) {
			Log.print( "Automatic depth correction to fault top @ 10 km" );
			eq.fault[i].depth = eq.fault[i].width/2 * sindeg(eq.fault[i].dip) + 10.e3;
		  }
		  else {
			Err.enable();
			return ierr;
		  }
		  ierr = eq.finalizeInput();
		}
		Err.enable();

    } else {

		// check fault parameters
		Err.disable();
		ierr = eq.finalizeInput();
		if( ierr ) {
		  i = ierr/10;
		  ierr = ierr - 10*i;
		  if( ierr != FLT_ERR_STRIKE ) {
			Err.enable();
			ierr = eq.finalizeInput();
			return ierr;
		  }
		  Log.print( "No strike on input: Employing effective symmetric source model" );
		  Err.enable();
		  if( eq.nfault > 1 ) return Err.post("symmetric source assumes only 1 fault");
		  eq.fault[0].strike = 0.;
		  effSymSource = 1;
		  ierr = eq.finalizeInput(); if(ierr) return ierr;
		}
		Err.enable();

    }

    // calculate uplift on a rectangular grid
    // set grid resolution, grid dimensions will be set automatically
    uZ.dx = DLon; uZ.dy = DLat;
    ierr = eq.calculate( uZ ); if(ierr) return ierr;
    
    if( effSymSource ) {
      // integrate for tsunami energy
      energy = 0.;
      for( j=0; j<uZ.ny; j++ ) {
        factLat = Dx*cosdeg(uZ.getY(0,j))*Dy;
        for( i=0; i<uZ.nx; i++ )
          energy += pow(uZ(i,j),2.)*factLat;
      }
      energy *= (1000*9.81/2);
      effRad = eq.fault[0].length/sqrt(2*M_PI);
      effMax = 1./effRad / sqrt(M_PI/2) / sqrt(1000*9.81/2) * sqrt(energy);
      Log.print( "Effective source radius: %g km,  max height: %g m", effRad/1000, effMax );
  
      // transfer uplift onto tsunami grid and define deformed area for acceleration
      for( i=0; i<uZ.nx; i++ ) {
        for( j=0; j<uZ.ny; j++ ) {
          dist = GeoDistOnSphere( uZ.getX(i,j),uZ.getY(i,j), eq.fault[0].lon,eq.fault[0].lat ) * 1000;
          if( dist < effRad ) 
            uZ(i,j) = effMax*cos(M_PI/2*dist/effRad);
          else
            uZ(i,j) = 0.;
        }
      }
      
    } // effective source
  
  } // src_type == fault

  // remove noise in the source
  absuzmax = uZ.getMaxAbsVal();

  if( (Par.ssh0ThresholdRel + Par.ssh0ThresholdAbs) != 0 ) {

    absuzmin = RealMax;
    if( Par.ssh0ThresholdRel != 0 ) absuzmin = Par.ssh0ThresholdRel*absuzmax;
    if( Par.ssh0ThresholdAbs != 0 && Par.ssh0ThresholdAbs < absuzmin ) absuzmin = Par.ssh0ThresholdAbs;

    for( i=0; i<uZ.nx; i++ ) {
      for( j=0; j<uZ.ny; j++ ) {
        if( fabs(uZ(i,j)) < absuzmin ) uZ(i,j) = 0;
      }
    }

  }

  // calculated (if needed) arrival threshold (negative value means it is relative)
  if( Par.sshArrivalThreshold < 0 ) Par.sshArrivalThreshold = absuzmax * fabs(Par.sshArrivalThreshold);

  // transfer uplift onto tsunami grid and define deformed area for acceleration
  Imin = NLon; Imax = 1; Jmin = NLat; Jmax = 1;
  /* FIXME: change loops */
  for( i=1; i<=NLon; i++ ) {
    for( j=1; j<=NLat; j++ ) {

      lon = getLon(i);
      lat = getLat(j);

      if( Node(idx(j,i), iD) != 0. )
        dz = Node(idx(j,i), iH) = uZ.getVal( lon,lat );
      else
        dz = Node(idx(j,i), iH) = 0.;

      if( fabs(dz) > Par.sshClipThreshold ) {
        Imin = My_min( Imin, i );
        Imax = My_max( Imax, i );
        Jmin = My_min( Jmin, j );
        Jmax = My_max( Jmax, j );
      }

    }
  }

  if( Imin == NLon ) return Err.post( "Zero initial displacement" );

  Imin = My_max( Imin - 2, 2 );
  Imax = My_min( Imax + 2, NLon-1 );
  Jmin = My_max( Jmin - 2, 2 );
  Jmax = My_min( Jmax + 2, NLat-1 );

  return 0;
}