//-----------------------------------------------------------------------
//
//            (c) Copyright 2006 by GATS, Inc.,
//     11864 Canon Blvd, Suite 101, Newport News VA 23606
//
//  All Rights Reserved. No part of this software or publication may be
//  reproduced, stored in a retrieval system, or transmitted, in any form
//  or by any means, electronic, mechanical, photocopying, recording, or
//  otherwise without the prior written permission of GATS, Inc.
//
//-----------------------------------------------------------------------
//
// Module:	RefracSimulation.cp
//
// Author:	Lance Deaver <l.e.deaver@gats-inc.com>
//
// Date:  	October 2, 2007
//
//-----------------------------------------------------------------------
//
//
//
//-----------------------------------------------------------------------
//
//-----------------------------------------------------------------------
// Include Files:
//-----------------------------------------------------------------------
//

#include "ConfigFile.h"
#include "Event.h"
#include "RefracSimulation.h"
#include "F77Radtran.h"
#include "F77Functions.h"
//#include "FORTRANarray.hpp"
#include "GATS_Utilities.hpp"

#include <cassert>
//
//-----------------------------------------------------------------------
// Defines and Macros:
//-----------------------------------------------------------------------
//
//#define THIS_MODULE	"RefracSimulation"


//
//-----------------------------------------------------------------------
// Global Variables:
//-----------------------------------------------------------------------
//

// there are several global routines and structures defined in FORTRAN code


//
//-----------------------------------------------------------------------
// Utility Routines:
//-----------------------------------------------------------------------
//

using namespace GATS_Utilities;

int RefracSimulation(Event& L0, Event& L1, Event& Tmp, Event& SD, ConfigFile& cf)
{


		
	static int mcel;
	static int ncel;
	static int mpout;

	static int c_one = 1;
	static int debug = 0;   		// 0 is a false (this is the debug flag)
	static int irf =  1;			// 1 is true always do refraction
	static int iplat = 0;  		// observer in space
	static int ic1 = -1;
	static int ic2 = 1;
	static  bool firstCall = true;
	static float wave;

	static std::vector<float> za, zt, pt, tt, gr, pout;
	static std::vector<double> zdub;
	static std::string  AltEventVarName,
		RefracEventVarName, sectionName ;

//	float lat83;
	float  re, tearth;
	int mr, nr;
	std::vector<double> refractionAng;
	std::vector<float> zr,pr,tr,tgr;



//  executable code starts here  
// initialization stuff done once
	if( firstCall  ) {
		sectionName="RefracSimulation";

		assert(cf.ValidEntry(sectionName , "WaveNumber") );
		assert(cf.ValidEntry(sectionName , "ImpactZ") );
		assert(cf.ValidEntry(sectionName, "AltEventVarName") );
		assert(cf.ValidEntry(sectionName, "RefracEventVarName") );

		wave =  cf.GetReal(sectionName , "WaveNumber" ) ;		
		za = cf.GetRealValues(sectionName,"ImpactZ");
		assert(monotonic_if(za.begin(), za.end(), std::greater<float>() ) );
 		zdub = std::vector<double>(za.begin(), za.end() );
                AltEventVarName = cf.GetStr(sectionName, "AltEventVarName");
                RefracEventVarName = cf.GetStr(sectionName, "RefracEventVarName");

		ncel=mcel=(int)za.size() ;
    		int it2 = mcel * mcel;
    		int it1 = it2 + mcel;
    		mpout = mcel * 11 + 3 + it1 * 12 + it2 * 3;
		pout.resize(mpout);
		zt.resize(mcel);
		pt.resize(mcel);
		tt.resize(mcel);
		gr.resize(2*mcel);	
// Call  KP_inita  so we can PROCEED
	
		::kp_inita__(&debug, &mcel, (int*)&pout[0]  )  ;

		firstCall = false;
	}


// Pull data from the L1 Event object.
  	EventVar<double> Z, P, T,  radCurvature ;
	L1.getEventVar("Z_Merged", Z);  //  merged NCEP and MSIS z,p,t profile in hydrostatics
	L1.getEventVar("P_Merged", P);
	L1.getEventVar("T_Merged", T);
	L1.getEventVar("L1_CurvatureRadius", radCurvature);
	re = (float)radCurvature[0];
		
// make floats from doubles
	zr = std::vector<float>(&Z[0], &Z[0]+Z.size() );
	pr = std::vector<float>(&P[0], &P[0]+P.size() );
	tr = std::vector<float>(&T[0], &T[0]+T.size() );	
	nr = mr=(int)zr.size();
// no temperature Gradients	
	tgr = std::vector<float>(2*mr, 0.);

// call path for the ray geometry

	::kp_mpath__( &mr,  &mcel,  &mpout,  &irf,  &c_one, &ncel, &c_one, &re, &ncel, &nr,
		&wave, &iplat, &za[0], &zr[0], &pr[0], &tr[0], &tgr[0],
 		&zt[0], &pt[0], &tt[0], &gr[0],&tearth, &pout[0] );

	
// Calculate the total Refraction angle along the limb path for each for the ray. (radians)
	refractionAng.clear();
	for(int ir=1; ir<= ncel; ++ir ) {
		refractionAng.push_back( ::refraction_( &ir, &ic1, &ic2, &pout[0] ) );	
	}
// make the EventVar's for this Band
	
	
	EventVar<double>  ev2( AltEventVarName , 
		sectionName+" Impact Altitudes", "km",
		&zdub[0], zdub.size() );

	L1.addEventVar( ev2 );

	EventVar<double> ev3(RefracEventVarName,
		sectionName+" Calculated Refraction Angle for "
		+ConvertToString(wave)+ " wavenumbers","radians",
		&refractionAng[0], refractionAng.size() );
	L1.addEventVar(ev3);
  	return 0;
}