/*
---- Hector the simulator ----
   A fast simulator of particles through generic beamlines.
   J. de Favereau, X. Rouby ~~~ hector_devel@cp3.phys.ucl.ac.be

        http://www.fynu.ucl.ac.be/hector.html

   Centre de Physique des Particules et de Phénoménologie (CP3)
   Université Catholique de Louvain (UCL)
*/

/// \file H_BeamLine.cc
/// \brief Reads external files and retrieves features of the real beam optical elements

// c++ #includes
#include <iostream>
#include <fstream>
#include <sstream>

// local #includes
#include "H_BeamLine.h"
#include "H_BeamLineParser.h"
#include "H_RectEllipticAperture.h"
#include "H_RectangularAperture.h"
#include "H_EllipticAperture.h"
#include "H_CircularAperture.h"
#include "H_RectangularDipole.h"
#include "H_SectorDipole.h"
#include "H_HorizontalQuadrupole.h"
#include "H_VerticalQuadrupole.h"
#include "H_HorizontalKicker.h"
#include "H_VerticalKicker.h"
#include "H_RectangularCollimator.h"
#include "H_Marker.h"

using namespace std;

// Caution : mad conventions for vertically(horizontally) focusing quadrupoles
// are opposite to Wille's. See fill() for more details.

H_BeamLine::H_BeamLine(const int si, const float length) : H_AbstractBeamLine(length){
	direction = (si >= abs(si)) ? 1 : -1;
	ips=0;
	ipx=0;
	ipy=0;
	iptx=0;
	ipty=0;
}

H_BeamLine::H_BeamLine(const H_BeamLine& beam) : H_AbstractBeamLine(beam) {
	direction = beam.direction;
	ips = beam.ips;
	ipx = beam.ipx;
	ipy = beam.ipy;
	iptx = beam.iptx;
	ipty = beam.ipty;
}

H_BeamLine& H_BeamLine::operator=(const H_BeamLine& beam) {
	if(this==&beam) return *this;
	direction = beam.direction;
	ips = beam.ips;
	ipx = beam.ipx;
    ipy = beam.ipy;
    iptx = beam.iptx;
    ipty = beam.ipty;
	return *this;
}

void H_BeamLine::findIP(const string filename) {
	findIP(filename,"IP5");
	return;
}

void H_BeamLine::findIP(const string filename, const string ipname) {
	// searches for the IP position in the extended table.
	ifstream tabfile(filename.c_str());
		if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl;
	bool found = false;
	int N_col=0;
	string headers[40];  // table of all column headers
	int header_type[40]; // table of all column types

	string temp_string;
	H_BeamLineParser e;
	istringstream curstring;

	while (	getline(tabfile,temp_string)) {
		curstring.clear(); // needed when using istringstream::str(string) several times !
		curstring.str(temp_string);

		// gets the features of each element
		if (found) {
			for (int col=0; col < N_col; col++) e.setProperties(curstring,header_type[col]);
			if(strstr(e.name.c_str(),ipname.c_str())) {
				ips = e.s; //e.printProperties();
				ipx = e.x;
				ipy = e.y;
				iptx = e.px*URAD;
				ipty = e.py*URAD;
			}
		} else if (strstr(temp_string.c_str(),"K0L")) {  //searches for the beginning of the element list.
			found = true;

			// reads the title of each column
			while(curstring.good()) { curstring >> headers[N_col]; if(headers[N_col]!="*") N_col++;}
			if(VERBOSE) cout << N_col << " columns identified" << endl;

			// identifies each column
			for (int col=0; col< N_col; col++)	header_type[col] =  column_identification(headers[col]);
			
			getline(tabfile,temp_string); // skip one line
		} // if(... "K0L")
	} // while (!eof)

	if (!found) cout << "\t ERROR ! IP not found." << endl;
	tabfile.close();
}

double* H_BeamLine::getIPProperties() {
	double* temp = new double[4];
	temp[0] = ipx;
	temp[1] = ipy;
	temp[2] = iptx;
	temp[3] = ipty;
	return temp;
}

void H_BeamLine::fill(const string filename) {
	fill(filename,1,"IP5");
	return;
}


void H_BeamLine::fill(const string filename, const int dir, const string ipname) {
	string headers[40];  // table of all column headers
	int header_type[40]; // table of all column types
	findIP(filename,ipname);
	ifstream tabfile(filename.c_str());
		if (! tabfile.is_open()) cout << "\t ERROR: I Can't open \"" << filename << "\"" << endl;
	if(VERBOSE) {
		cout<<"Using file : "<< filename <<" in the "<<((direction>0)?"positive":"negative")<<" direction."<<endl;
		cout<<"IP was found at position : "<<ips<<endl<<endl;
	}
	bool found = false; // found = true when the beginning of the element list is reached
	int type, N_col=0;
	float previous_betax =0, previous_betay=0; // needed when direction ==1
	float previous_dx =0, previous_dy=0; // needed when direction ==1
	float previous_x =0, previous_y=0; // needed when direction ==1

	string temp_string;
	H_BeamLineParser e;
	istringstream curstring;
	H_OpticalElement *el = 0;
	H_Aperture * ap = 0;


	while (getline(tabfile,temp_string)) {
		curstring.clear(); // needed when using several tims istringstream::str(string)
		curstring.str(temp_string);

        // gets the features of each element
		if (found) {
			// reads each column
			for (int col=0; col < N_col; col++) {
				e.setProperties(curstring,header_type[col]);
			}
			//e.printProperties();
			
			type =0; //init
			if(e.k0l!=0) { if(strstr(e.name.c_str(),"MB.")) type = SDIPOLE; else type = RDIPOLE;} //all SBEND seems to be called MB.xxxxx.xx
			else if(e.k1l!=0) type = (e.k1l>0)?HQUADRUPOLE:VQUADRUPOLE;
			else if(e.hkick!=0) type = HKICKER;
			else if(e.vkick!=0) type = VKICKER;
			else if(strstr(e.name.c_str(),"DRIFT")) type=DRIFT;
			else if((strstr(e.name.c_str(),"DFB"))||(strstr(e.name.c_str(),"TAS"))||(strstr(e.name.c_str(),"TAN"))||(strstr(e.name.c_str(),"TCL"))) type = RCOLLIMATOR;
			else if(strstr(e.name.c_str(),ipname.c_str())) { type=IP; }
			else type=MARKER;
			
			e.s = (e.s > ips) ? (e.s -ips - e.l)*(direction) : (e.s-ips)*(direction);
			
			el=0; //init
			switch(type) {
				case RDIPOLE:     { el = new H_RectangularDipole(e.name.c_str(),e.s,dir*e.k0l/e.l,e.l); } break;
				case SDIPOLE:     { el = new H_SectorDipole(e.name.c_str(),e.s,dir*e.k0l/e.l,e.l); } break;
				case VQUADRUPOLE: { el = new H_VerticalQuadrupole(e.name.c_str(),e.s,-(e.k1l/e.l),e.l); } break;
				case HQUADRUPOLE: { el = new H_HorizontalQuadrupole(e.name.c_str(),e.s,-(e.k1l/e.l),e.l); } break;
				case VKICKER:     { el = new H_VerticalKicker(e.name.c_str(),e.s,e.vkick,e.l); } break;
				case HKICKER:     { el = new H_HorizontalKicker(e.name.c_str(),e.s,e.hkick,e.l); } break;
				case RCOLLIMATOR: { el = new H_RectangularCollimator(e.name.c_str(),e.s,e.l); } break;
				case DRIFT: {previous_betax = e.betx; previous_betay = e.bety; previous_dx = e.dx; previous_dy = e.dy; previous_x = e.x; previous_y = e.y;} break;
				case IP: { el = new H_Marker(e.name.c_str(),e.s); } break;
				default: break;
			}

			ap = 0; //init
			if(e.aper_1 !=0 || e.aper_2 !=0 || e.aper_3 !=0 || e.aper_4 != 0) {
				e.aper_1 *= URAD;
				e.aper_2 *= URAD;
				e.aper_3 *= URAD;
				e.aper_4 *= URAD; // in [m] in the tables !

				if(strstr(e.apertype.c_str(),"RECTELLIPSE")) ap = new H_RectEllipticAperture(e.aper_1,e.aper_2,e.aper_3,e.aper_4,0,0);
				else if(strstr(e.apertype.c_str(),"CIRCLE")) ap = new H_CircularAperture(e.aper_1,0,0);
				else if(strstr(e.apertype.c_str(),"RECTANGLE")) ap = new H_RectangularAperture(e.aper_1,e.aper_2,0,0);
				else if(strstr(e.apertype.c_str(),"ELLIPSE")) ap = new H_EllipticAperture(e.aper_1,e.aper_2,0,0);
			}

			if(el!=0) {
				if (direction<0) {
					el->setBetaX(e.betx);
					el->setBetaY(e.bety);
					el->setDX(e.dx);
					el->setDY(e.dy);
					el->setRelX(e.x);
					el->setRelY(e.y);
				} else {
					el->setBetaX(previous_betax);
					el->setBetaY(previous_betay);
					el->setDX(previous_dx);
					el->setDY(previous_dy);
					el->setRelX(previous_x);
					el->setRelY(previous_y);
				}
				if(ap!=0) { 
					el->setAperture(ap); 
				//		delete ap; // ap deleted in H_AbstractBeamLine::~H_AbstractBeamLine
				}
	
				/// Parses all the elements, but only keeps the ones from the IP till the desired length
				if(e.s>=0 && e.s<beam_length) add(*el);

			// delete el; // el deleted in H_AbstractBeamLine::~H_AbstractBeamLine
			} // if(el!=0)
		} // if (found)
		else if(strstr(temp_string.c_str(),"K0L")) { // if (!found)
		//searches for the beginning of the element list.
			found = true;

			// reads the title of each column
			while(curstring.good()) { curstring >> headers[N_col]; if(headers[N_col]!="*") N_col++;}
			if(VERBOSE) cout << N_col << " columns identified" << endl;

			// identifies each column
			for (int col=0; col< N_col; col++) {
				header_type[col] =  column_identification(headers[col]);
			}
			getline(tabfile,temp_string);
		} // if temp_string <-> "K0L"

	} // while (!eof)
	tabfile.close();
}