// $Id: rtfromfpar.cpp,v 1.1 2009/02/09 20:57:21 thernis Exp $


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iomanip>
#include <string>
#include <sstream>
#include <vector>

#if (FORCE_CCFITS!=1)  // Don't use CCfits if not explicitly requested when configure
#undef HAVE_LIBCCFITS
#endif

#ifdef HAVE_LIBCCFITS
#warning "Using CCfits"
#include <CCfits/CCfits>
#endif
#include "Cvec.h"
#include "Cmat.h"
#include "Cbasis.h"
#include "CModelBase.h"
#include "raytrace.h"

using namespace std;

int main(int argc, char **argv);
#ifdef HAVE_LIBCCFITS
void saveResultsInFits(rtparam fp,int nbtotpix);
void fillinheader(CCfits::FITS* pFits,rtparam fp);
void transferinvalarray(float* pim,int nbtotpix,valarray<float> &array);
#endif
//int savedatfile(char* filename,float *pim,int nbtotpix);
//int savefitsfile(char* filename,float *pim,int nbtotpix,long nbaxis,long* naxes);
//void printerror(int status);

//! Get the parameters from a file and do the raytracing
int main(int argc, char **argv)
{
	// -- get the filename of the parameter file
	argv++;
	ostringstream ostrFilename;
	ostrFilename << *argv;
		
	// -- read the parameter file
	cout << "opening file : " << ostrFilename.str() << endl;
	ifstream myFile(ostrFilename.str().c_str(),ios::in);
	
	rtparam fp;
	char junk[100];
	
	// -- the file should always be formated the same
	// -- modelid
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int modelid=atoi(junk);
	fp.pmodelid=&modelid;
	// -- imsize
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int is,js;
	is=atoi(junk);
	fp.pis=&is;
	myFile.getline(junk,100);
	cout << junk << endl;
	js=atoi(junk);
	fp.pjs=&js;
	// -- fovpix
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float fovpix=atof(junk);
	fp.pfovpix=&fovpix;
	// -- obspos
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float obspos[3];
	obspos[0]=atof(junk);
	fp.pobspos=&obspos[0];
	myFile.getline(junk,100);
	cout << junk << endl;
	obspos[1]=atof(junk);
	myFile.getline(junk,100);
	cout << junk << endl;
	obspos[2]=atof(junk);
	// -- obsang
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float obsang[3];
	obsang[0]=atof(junk);
	fp.pobsang=&obsang[0];
	myFile.getline(junk,100);
	cout << junk << endl;
	obsang[1]=atof(junk);
	myFile.getline(junk,100);
	cout << junk << endl;
	obsang[2]=atof(junk);
	// -- nepos
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float nepos[3];
	nepos[0]=atof(junk);
	fp.pnepos=&nepos[0];
	myFile.getline(junk,100);
	cout << junk << endl;
	nepos[1]=atof(junk);
	myFile.getline(junk,100);
	cout << junk << endl;
	nepos[2]=atof(junk);
	// -- neang
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float neang[3];
	neang[0]=atof(junk);
	fp.pneang=&neang[0];
	myFile.getline(junk,100);
	cout << junk << endl;
	neang[1]=atof(junk);
	myFile.getline(junk,100);
	cout << junk << endl;
	neang[2]=atof(junk);
	// -- losnbp
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int losnbp=atoi(junk);
	fp.plosnbp=&losnbp;
	// -- losrange
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	float losrange[2];
	losrange[0]=atof(junk);
	fp.plosrange=&losrange[0];
	myFile.getline(junk,100);
	cout << junk << endl;
	losrange[1]=atof(junk);
	// -- modparam
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	fp.pmparam=new float;
	*(fp.pmparam)=atof(junk);
	// -- pofinteg
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int pofinteg=atoi(junk);
	fp.ppofinteg=&pofinteg;
	// -- quiet
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int quiet=atoi(junk);
	fp.pquiet=&quiet;
	// -- neonly
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	int neonly=atoi(junk);
	fp.pneonly=&neonly;
	// -- neonly
	myFile.getline(junk,100);
	cout << junk << endl;
	myFile.getline(junk,100);
	cout << junk << endl;
	fp.Hlonlat=new float[3];
	fp.Hlonlat[0]=atof(junk);
	myFile.getline(junk,100);
	cout << junk << endl;
	fp.Hlonlat[1]=atof(junk);
	cout << junk << endl;
	fp.Hlonlat[2]=atof(junk);

	myFile.close();

	// ---- init multiple output from rt
	int nbtotpix=(*fp.pis)*(*fp.pjs);
	// -- btot
	fp.pbtot= new float[nbtotpix];
	fp.pbpol= new float[nbtotpix];
	fp.pnetot= new float[nbtotpix];
	fp.prhoim= new float[nbtotpix];
	fp.pmmlon= new float[2];
	fp.pmmlat= new float[2];
	fp.prrr= new float[nbtotpix];
	fp.pcrpix=new float[2];
	fp.proi=new int[nbtotpix];

	// -- set the points of the ROI to 1
	for(int i=0;i<nbtotpix;i++) {
		*(fp.proi + i)=1;
	}

	cout << "I'm here ! " << endl;
	
	cout << "fp.pis : " << *fp.pis << endl;
	cout << "fp.pjs : " << *fp.pjs << endl;
	cout << "fp.ppofinteg : " << *(fp.ppofinteg) << endl;
	cout << "fp.plosrange[0] : " << fp.plosrange[0] << endl;
	cout << "fp.plosrange[1] : " << fp.plosrange[1] << endl;
	cout << "fp.plosnbp : " << *(fp.plosnbp) << endl;
	cout << "fp.pobspos[0] : " << fp.pobspos[0] << endl;
	cout << "fp.pobspos[1] : " << fp.pobspos[1] << endl;
	cout << "fp.pobspos[2] : " << fp.pobspos[2] << endl;
	
	// -- call the raytracing routine
	int rtxitstat=raytracemain(fp);
	
	cout << "Exit status raytracemain : " << rtxitstat << endl;
	
	// ---- save in fits format
#ifdef HAVE_LIBCCFITS
	saveResultsInFits(fp,nbtotpix);
#else
 cout << "Cannot save the image as a fits file: Raytrace not built with CCFITS." << endl;
#endif
	// -- clean up the mess
	delete [] fp.pmparam;
	delete [] fp.pbtot;
	delete [] fp.pbpol;
	delete [] fp.pnetot;
	delete [] fp.prhoim;
	delete [] fp.pmmlon;
	delete [] fp.pmmlat;
	delete [] fp.prrr;
	delete [] fp.pcrpix;
	delete [] fp.proi;
	delete [] fp.Hlonlat;
	cout << "I'm there ! " << endl;

	return 1;
}


//! Save the raytrace results in fits format
#ifdef HAVE_LIBCCFITS
void saveResultsInFits(rtparam fp,int nbtotpix)
{
	// -- define axis number and size
	long nbaxis=2;
	long naxes[2];
	naxes[0]=*fp.pis;naxes[1]=*fp.pjs;
	vector<long> Vnaxes;
	Vnaxes.push_back(naxes[0]);
	Vnaxes.push_back(naxes[1]);
	long fpixel(1);
	valarray<float> array(nbtotpix);

	
	// --- save btot
	const string btotfn("!btot.fts");
	// -- open file
	auto_ptr<CCfits::FITS> pFits(0);
	pFits.reset(new CCfits::FITS(btotfn,FLOAT_IMG,nbaxis,naxes));
	// -- write the image
	transferinvalarray(fp.pbtot,nbtotpix,array);
	string headername("HEADER");
	CCfits::ExtHDU* imageExt=pFits->addImage(headername,FLOAT_IMG,Vnaxes);
	imageExt->write(fpixel,nbtotpix,array);
	// -- build the header
	fillinheader(pFits.get(),fp);
	// -- save the header
	pFits->pHDU().write(fpixel,nbtotpix,array);
	pFits->pHDU().addKey("SIMUTYPE","Total Brightness","");
	// -- close the file
	pFits->destroy();

	// --- save bpol
	const string bpolfn("!bpol.fts");
	// -- open file
	pFits.reset(new CCfits::FITS(bpolfn,FLOAT_IMG,nbaxis,naxes));
	// -- write the image
	transferinvalarray(fp.pbpol,nbtotpix,array);
	imageExt=pFits->addImage(headername,FLOAT_IMG,Vnaxes);
	imageExt->write(fpixel,nbtotpix,array);
	// -- build the header
	fillinheader(pFits.get(),fp);
	pFits->pHDU().addKey("SIMUTYPE","Polarized Brightness","");
	// -- save the header
	pFits->pHDU().write(fpixel,nbtotpix,array);
	// -- close the file
	pFits->destroy();

	// --- save netot
	const string netotfn("!netot.fts");
	// -- open file
	pFits.reset(new CCfits::FITS(netotfn,FLOAT_IMG,nbaxis,naxes));
	// -- write the image
	transferinvalarray(fp.pnetot,nbtotpix,array);
	imageExt=pFits->addImage(headername,FLOAT_IMG,Vnaxes);
	imageExt->write(fpixel,nbtotpix,array);
	// -- build the header
	fillinheader(pFits.get(),fp);
	pFits->pHDU().addKey("SIMUTYPE","Integrated Ne","");
	// -- save the header
	pFits->pHDU().write(fpixel,nbtotpix,array);
	// -- close the file
	pFits->destroy();	
	
	
	
}


//! Fill in the fits header with the raytrace parameters
void fillinheader(CCfits::FITS* pFits,rtparam fp)
{
	pFits->pHDU().addKey("FOVPIX",*fp.pfovpix,"rad");
	pFits->pHDU().addKey("OBSPOSX",fp.pobspos[0],"Rsun");
	pFits->pHDU().addKey("OBSPOSY",fp.pobspos[1],"Rsun");
	pFits->pHDU().addKey("OBSPOSZ",fp.pobspos[2],"Rsun");
	pFits->pHDU().addKey("OBSANGX",fp.pobsang[0],"rad");
	pFits->pHDU().addKey("OBSANGY",fp.pobsang[1],"rad");
	pFits->pHDU().addKey("OBSANGZ",fp.pobsang[2],"rad");
	pFits->pHDU().addKey("NEPOSX",fp.pnepos[0],"Rsun");
	pFits->pHDU().addKey("NEPOSY",fp.pnepos[1],"Rsun");
	pFits->pHDU().addKey("NEPOSZ",fp.pnepos[2],"Rsun");
	pFits->pHDU().addKey("NEANGX",fp.pneang[0],"rad");
	pFits->pHDU().addKey("NEANGY",fp.pneang[1],"rad");
	pFits->pHDU().addKey("NEANGZ",fp.pneang[2],"rad");
	pFits->pHDU().addKey("LOSNBP",*fp.plosnbp,"points");
	pFits->pHDU().addKey("LOSRNG1",fp.plosrange[0],"Rsun");
	pFits->pHDU().addKey("LOSRNG2",fp.plosrange[1],"Rsun");
	pFits->pHDU().addKey("MODELID",*fp.pmodelid,"");
}


//! transfer the image in a valarray
void transferinvalarray(float* pim,int nbtotpix,valarray<float> &array)
{
    // -- put the image in the array
	float *pscan;
	pscan=pim;
	for (int i=0;i<nbtotpix;i++) array[i]=*(pscan++);
}


#endif


/*
* $Log: rtfromfpar.cpp,v $
* Revision 1.1  2009/02/09 20:57:21  thernis
* Put some old code here in case I need that one day...
*
* Revision 1.4  2007/05/14 17:19:41  thernis
* Add CVS id and log in all files
*
*/