#include <stdlib.h>

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>
#include <complex>

#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>

using namespace std;
using namespace boost;

typedef complex<double> Complex;

class terminate_program {
	private:
		string msg;
	public:
		terminate_program() { msg=string(); }
		terminate_program(const char *e) { msg=string(e); }
		terminate_program(const string e) { msg=e; }
		string read() { return msg; }
		bool ifmsg() { if(msg.size()) return true; else return false; }
		terminate_program &operator+=(const string e) {
			this->msg=this->msg+e;
			return *this;
		}
};

class udm{
	private:
		double Nvc, Eg, Eh;        // parameters of interband transitions (eV^2, eV, eV)
		bool he;                   // check if higher energy itransitions are included in model
		double Nvx, Ex;            // parameters of higher energy itransitions (eV^2, eV)
		bool ut;                   // check if Urbach (exponential) tail is included in model
		double Nut, Eu;            // parameters of Urbach (exponential) tail (eV^2, eV)
		unsigned nex;              // number of excitons
		bool gex;                  // Gaussian extension of excitons
		double A0;                 // amplitude of symmetrical interband transitions term (usually 1, only with excitons)
		double *Aex;               // amplitude of excitonic terms of interband transitions
		double *Eex;               // mean energy of excitonic terms (eV)
		double *Bex;               // broadening of excitonic terms (eV)
		double *Gex;               // Gaussian part of exitonic terms in range [0,1]
		bool dEg;                  // doubled bang gap;
		double Ep;                 // mean phonon energy doubling gap
		unsigned nloc;             // number of peak corresponding transitions on localized states
		double *Nloc;              // transition strength of localized state transition peks (eV^2)
		double *Eloc;              // mean energy of localized state transition peks (eV)
		double *Bloc;              // broadening of localized state transition peks (eV)
		unsigned nph;              // number of phonon peaks
		double *Nph;               // transition strengt of phonon peaks (eV^2)
		double *nuph;              // mean energy (wavenumber) of phonon peaks (1/cm)
		double *betaph;            // FWHM of phonon peaks (1/cm)
		string notices;            // text notices
	public:
		udm(string parfile);       // constructor reading model specification and parameter values from <parfile>
		void arrange(void);        // arrange of numbered parameters
		void print(void);          // print parfile to stdout
		Complex df(double E);      // calculation of dielectric function for photon energy E (eV)
};

udm::udm(string parfile)
{
	ifstream in(parfile.c_str());
	if(in.fail()) {
		terminate_program err("can't open file ");
		err+=parfile;
		throw err;
	}

	vector<string> pars;
	pars.push_back("Nvc");
	pars.push_back("Eg");
	pars.push_back("Eh");
	he=false;
	ut=false;
	gex=false;
	dEg=false;
	nex=0;
	nloc=0;
	nph=0;

	string attributes;
	getline(in,attributes);

	vector<string> vatt;
	split( vatt, attributes, is_any_of( ":" ) );
	for (vector<string>::iterator att = vatt.begin(); att != vatt.end(); att++) {
		if( *att=="" ) ; else
		if( *att=="dEg" ) {
			pars.push_back("Ep");
			dEg=true;
		} else
		if( *att=="he" ) {
			pars.push_back("Nvx");
			pars.push_back("Ex");
			he=true;
		} else
		if( *att=="ut" ) {
			pars.push_back("Nut");
			pars.push_back("Eu");
			ut=true;
		} else
		if( att->substr(0,2)=="ex" || att->substr(0,3)=="gex" ) {
			if( att->substr(0,3)=="gex" ) gex=true;
			size_t p = att->find('=');
			string nexs=att->substr(p+1);
			trim(nexs);
			try{ nex=lexical_cast<unsigned>(nexs);}
			catch(bad_lexical_cast) {
				terminate_program err="attributes "+attributes+": can't covert \""+nexs+"\" to number";
				throw err;
			}
			if(nex) pars.push_back("A0");
			for(unsigned i=1;i<=nex;i++) {
				pars.push_back("Aex"+lexical_cast<string>(i));
				pars.push_back("Eex"+lexical_cast<string>(i));
				pars.push_back("Bex"+lexical_cast<string>(i));
				if(gex) pars.push_back("Gex"+lexical_cast<string>(i));
			}
			Aex = new double[nex];
			Eex = new double[nex];
			Bex = new double[nex];
			if(gex) Gex = new double[nex];
		} else
		if( att->substr(0,3)=="loc" ) {
			size_t p = att->find('=');
			string nlocs=att->substr(p+1);
			trim(nlocs);
			try{ nloc=lexical_cast<unsigned>(nlocs);}
			catch(bad_lexical_cast) {
				terminate_program err="attributes "+attributes+": can't covert \""+nlocs+"\" to number";
				throw err;
			}
			for(unsigned i=1;i<=nloc;i++) {
				pars.push_back("Nloc"+lexical_cast<string>(i));
				pars.push_back("Eloc"+lexical_cast<string>(i));
				pars.push_back("Bloc"+lexical_cast<string>(i));
			}
			Nloc = new double[nloc];
			Eloc = new double[nloc];
			Bloc = new double[nloc];
		} else
		if( att->substr(0,2)=="ph" ) {
			size_t p = att->find('=');
			string nphs=att->substr(p+1);
			trim(nphs);
			try{ nph=lexical_cast<unsigned>(nphs);}
			catch(bad_lexical_cast) {
				terminate_program err="attributes "+attributes+": can't covert \""+nphs+"\" to number";
				throw err;
			}
			for(unsigned i=1;i<=nph;i++) {
				pars.push_back("Nph"+lexical_cast<string>(i));
				pars.push_back("nuph"+lexical_cast<string>(i));
				pars.push_back("betaph"+lexical_cast<string>(i));
			}
			Nph = new double[nph];
			nuph = new double[nph];
			betaph = new double[nph];
		} else {
			terminate_program err="attributes "+attributes+": wrong attribute \""+*att+"\" of Universal model";
			throw err;
		}
	}

	vector<bool> parb(pars.size(),false);

	for(unsigned i=0;i<pars.size();i++) {
		string line;
		getline(in,line);

		if(in.fail()) {
			terminate_program err="few number of parameters\ndefine parameters:";
			for(int j=0;j<pars.size();j++) if(!parb[j]) err+=" "+pars[j];
			throw err;
		}

		size_t p=line.find('=');
		if( p == string::npos ) {
			terminate_program err="wrong "+lexical_cast<string>(i+1)+".line ("+line+"): can't find '='";
			throw err;
		}
		string parname=line.substr(0,p);
		trim(parname);
		unsigned j=0;
		for(;j<pars.size();j++) if( pars[j]==parname ) break;
		if(j==pars.size()) {
			terminate_program err="unknown parameter "+parname;
			throw err;
		}
		if(parb[j]) {
			terminate_program err="wrong "+lexical_cast<string>(i+1)+".line: parameter "+parname+" is already defined";
			throw err;
		}
		parb[j]=true;
		string values=line.substr(p+1);
		trim(values);
		double value;
		try{value=lexical_cast<double>(values);}
		catch(bad_lexical_cast) {
			terminate_program err="wrong "+lexical_cast<string>(i+1)+".line ("+line+"): can't covert \""+values+"\" to number";
			throw err;
		}
		if(parname=="Nvc") {Nvc=value; continue;}
		if(parname=="Eg") {Eg=value; continue;}
		if(parname=="Eh") {Eh=value; continue;}
		if(dEg) {
			if(parname=="Ep") {Ep=value; continue;}
		}
		if(nex) {
			if(parname=="A0") {A0=value; continue;}
			for(j=1;j<=nex;j++) if(parname=="Aex"+lexical_cast<string>(j)) {Aex[j-1]=value; break;}
			if(j<=nex) continue;
			for(j=1;j<=nex;j++) if(parname=="Eex"+lexical_cast<string>(j)) {Eex[j-1]=value; break;}
			if(j<=nex) continue;
			for(j=1;j<=nex;j++) if(parname=="Bex"+lexical_cast<string>(j)) {Bex[j-1]=value; break;}
			if(j<=nex) continue;
			if(gex) for(j=1;j<=nex;j++) if(parname=="Gex"+lexical_cast<string>(j)) {Gex[j-1]=value; break;}
			if(j<=nex) continue;
		} else A0=1;
		if(he) {
			if(parname=="Nvx") {Nvx=value; continue;}
			if(parname=="Ex") {Ex=value; continue;}
		}
		if(ut) {
			if(parname=="Nut") {Nut=value; continue;}
			if(parname=="Eu") {Eu=value; continue;}
		}
		if(nloc) {
			for(j=1;j<=nloc;j++) if(parname=="Nloc"+lexical_cast<string>(j)) {Nloc[j-1]=value; break;}
			if(j<=nloc) continue;
			for(j=1;j<=nloc;j++) if(parname=="Eloc"+lexical_cast<string>(j)) {Eloc[j-1]=value; break;}
			if(j<=nloc) continue;
			for(j=1;j<=nloc;j++) if(parname=="Bloc"+lexical_cast<string>(j)) {Bloc[j-1]=value; break;}
			if(j<=nloc) continue;
		}
		if(nph) {
			for(j=1;j<=nph;j++) if(parname=="Nph"+lexical_cast<string>(j)) {Nph[j-1]=value; break;}
			if(j<=nph) continue;
			for(j=1;j<=nph;j++) if(parname=="nuph"+lexical_cast<string>(j)) {nuph[j-1]=value; break;}
			if(j<=nph) continue;
			for(j=1;j<=nph;j++) if(parname=="betaph"+lexical_cast<string>(j)) {betaph[j-1]=value; break;}
			if(j<=nph) continue;
		}

		terminate_program err="inconsistency in parsing of "+parname+" parameter";
		throw err;
	}
	
	do {
		string line;
		getline(in,line);
		
		if(in.fail()) return;
		notices += line+"\n";
	} while(true);
}

static unsigned sort_peaks(unsigned n,double *E,double *A,double *B,double *G)
{
	for(unsigned i=0;i<n;i++) if(A[i]==0.0) {
		n--;
		for(unsigned j=i;j<n;j++) {
			A[j]=A[j+1];
			E[j]=E[j+1];
			B[j]=B[j+1];
			if(G) G[j]=G[j+1];
		}
		i--;
	}
	
	bool sw=true;
	
	while(sw) {
		sw=false;
		for(int i=1;i<n;i++) if(E[i-1]>E[i]) {
			sw=true;
			double t=E[i];
			E[i]=E[i-1];
			E[i-1]=t;
			t=A[i];
			A[i]=A[i-1];
			A[i-1]=t;
			t=B[i];
			B[i]=B[i-1];
			B[i-1]=t;
			if(G) {
				t=G[i];
				G[i]=G[i-1];
				G[i-1]=t;
			}
		}
	}
	
	return n;
}

void udm::arrange(void)
{
	if(nex>1) {
		if(gex) nex=sort_peaks(nex,Eex,Aex,Bex,Gex);
		else nex=sort_peaks(nex,Eex,Aex,Bex,NULL);
	}
	if(nloc>1) nloc=sort_peaks(nloc,Eloc,Nloc,Bloc,NULL);
	if(nph>1) nph=sort_peaks(nph,nuph,Nph,betaph,NULL);
}

void udm::print(void)
{
	string separator;
	if(he) {
		cout << "he";
		separator=":";
	}
	if(nex) {
		cout << separator << "ex=" << nex;
		separator=":";
	}
	if(dEg) {
		cout << separator << "dEg";
		separator=":";
	}
	if(ut) {
		cout << separator << "ut";
		separator=":";
	}
	if(nloc) {
		cout << separator << "loc=" << nloc;
		separator=":";
	}
	if(nph) {
		cout << separator << "ph=" << nph;
	}
	cout << endl;

	cout << "Nvc = " << std::setprecision(16) << Nvc << endl;
	cout << "Eg = " << std::setprecision(16) << Eg << endl;
	cout << "Eh = " << std::setprecision(16) << Eh << endl;
	if(nex) {
		cout << "A0 = " << std::setprecision(16) << A0 << endl;
		for(unsigned i=0;i<nex;i++) {
			cout << "Aex" << i+1 << " = " << std::setprecision(16) << Aex[i] << endl;
			cout << "Eex" << i+1 << " = " << std::setprecision(16) << Eex[i] << endl;
			cout << "Bex" << i+1 << " = " << std::setprecision(16) << Bex[i] << endl;
			if(gex) cout << "Gex" << i+1 << " = " << std::setprecision(16) << Gex[i] << endl;
		}
	}
	if(he) {
		cout << "Nvx = " << std::setprecision(16) << Nvx << endl;
		cout << "Ex = " << std::setprecision(16) << Ex << endl;
	}
	if(ut) {
		cout << "Nut = " << std::setprecision(16) << Nut << endl;
		cout << "Eu = " << std::setprecision(16) << Eu << endl;
	}
	if(nloc) {
		for(unsigned i=0;i<nloc;i++) {
			cout << "Nloc" << i+1 << " = " << std::setprecision(16) << Nloc[i] << endl;
			cout << "Eloc" << i+1 << " = " << std::setprecision(16) << Eloc[i] << endl;
			cout << "Bloc" << i+1 << " = " << std::setprecision(16) << Bloc[i] << endl;
		}
	}
	if(nph) {
		for(unsigned i=0;i<nph;i++) {
			cout << "Nph" << i+1 << " = " << std::setprecision(16) << Nph[i] << endl;
			cout << "nuph" << i+1 << " = " << std::setprecision(16) << nuph[i] << endl;
			cout << "betaph" << i+1 << " = " << std::setprecision(16) << betaph[i] << endl;
		}
	}
	if(notices.size()) cout << notices << endl;
}


#define M_SQRTPI       1.77245385090551602729816748334      /* sqrt(pi) */
#define M_EULER        0.57721566490153286060651209008      /* Euler constant */
#define Planck         6.6260693e-34
#define light_speed    299792458
#define el_charge      1.60217653e-19
#define tnu_to_E       (1e2 * Planck * light_speed / el_charge)
#define nu_to_E        (1e12 * Planck / el_charge)
#define lambda_E       (1e9 * Planck * light_speed / el_charge)

/* The interval is always [-1..1], it should be precomputed to avoid doing the transform here. */
template <typename T, unsigned int N>
class ChebyshevSeries {
	public:
		const unsigned int order_long;
		const unsigned int order_double;
		const unsigned int order_single;
		const T c[N];

	private:
		inline T eval_terms(const T x, unsigned int nterms) const
		{
			T d  = 0.0;
			T dd = 0.0;
			T x2 = 2*x;

			for (unsigned int j = nterms; j >= 1; j--) {
				T temp = d;
				d = x2*d - dd + c[j];
				dd = temp;
			}

			return x*d - dd + 0.5*c[0];
		}

	public:
		inline T eval_long(const T x) const { return eval_terms(x, order_long); }
		inline T eval_double(const T x) const { return eval_terms(x, order_double); }
		inline T eval_single(const T x) const { return eval_terms(x, order_single); }
};
#define spi2 0.8862269254527580136490837416705725913990

/* Library libcerf:
 *   Compute complex error functions, based on a new implementation of
 *   Faddeeva's w_of_z. Also provide Dawson and Voigt functions.
 *
 * File im_w_of_x.c:
 *   Compute scaled Dawson integral im_w_of_x(x) = 2*dawson(x)/sqrt(pi),
 *   equivalent to the imaginary part of the Faddeeva function w(x) for real x.
 * 
 * Copyright:
 *   (C) 2012 Massachusetts Institute of Technology
 *   (C) 2013 Forschungszentrum Jülich GmbH
 * 
 * Licence:
 *   Permission is hereby granted, free of charge, to any person obtaining
 *   a copy of this software and associated documentation files (the
 *   "Software"), to deal in the Software without restriction, including
 *   without limitation the rights to use, copy, modify, merge, publish,
 *   distribute, sublicense, and/or sell copies of the Software, and to
 *   permit persons to whom the Software is furnished to do so, subject to
 *   the following conditions:
 * 
 *   The above copyright notice and this permission notice shall be
 *   included in all copies or substantial portions of the Software.
 * 
 *   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
 *   EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
 *   MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
 *   NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
 *   LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
 *   OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
 *   WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. 
 *
 * Authors:
 *   Steven G. Johnson, Massachusetts Institute of Technology, 2012, core author
 *   Joachim Wuttke, Forschungszentrum Jülich, 2013, package maintainer
 *
 * Website:
 *   http://apps.jcns.fz-juelich.de/libcerf
 *
 * Revision history:
 *   ../CHANGELOG
 *
 * Manual page:
 *   man 3 im_w_of_x
 */

/******************************************************************************/
/* Lookup-table for Chebyshev polynomials for smaller |x|                     */
/******************************************************************************/

	static double
w_im_y100(double y100, double x)
{

	// Steven G. Johnson, October 2012.

	// Given y100=100*y, where y = 1/(1+x) for x >= 0, compute w_im(x).

	// Uses a look-up table of 100 different Chebyshev polynomials
	// for y intervals [0,0.01], [0.01,0.02], ...., [0.99,1], generated
	// with the help of Maple and a little shell script.
	// This allows the Chebyshev polynomials to be of significantly lower
	// degree (about 1/30) compared to fitting the whole [0,1] interval
	// with a single polynomial.

	switch ((int) y100) {
		case 0: {
				double t = 2*y100 - 1;
				return 0.28351593328822191546e-2 + (0.28494783221378400759e-2 + (0.14427470563276734183e-4 + (0.10939723080231588129e-6 + (0.92474307943275042045e-9 + (0.89128907666450075245e-11 + 0.92974121935111111110e-13 * t) * t) * t) * t) * t) * t;
			}
		case 1: {
				double t = 2*y100 - 3;
				return 0.85927161243940350562e-2 + (0.29085312941641339862e-2 + (0.15106783707725582090e-4 + (0.11716709978531327367e-6 + (0.10197387816021040024e-8 + (0.10122678863073360769e-10 + 0.10917479678400000000e-12 * t) * t) * t) * t) * t) * t;
			}
		case 2: {
				double t = 2*y100 - 5;
				return 0.14471159831187703054e-1 + (0.29703978970263836210e-2 + (0.15835096760173030976e-4 + (0.12574803383199211596e-6 + (0.11278672159518415848e-8 + (0.11547462300333495797e-10 + 0.12894535335111111111e-12 * t) * t) * t) * t) * t) * t;
			}
		case 3: {
				double t = 2*y100 - 7;
				return 0.20476320420324610618e-1 + (0.30352843012898665856e-2 + (0.16617609387003727409e-4 + (0.13525429711163116103e-6 + (0.12515095552507169013e-8 + (0.13235687543603382345e-10 + 0.15326595042666666667e-12 * t) * t) * t) * t) * t) * t;
			}
		case 4: {
				double t = 2*y100 - 9;
				return 0.26614461952489004566e-1 + (0.31034189276234947088e-2 + (0.17460268109986214274e-4 + (0.14582130824485709573e-6 + (0.13935959083809746345e-8 + (0.15249438072998932900e-10 + 0.18344741882133333333e-12 * t) * t) * t) * t) * t) * t;
			}
		case 5: {
				double t = 2*y100 - 11;
				return 0.32892330248093586215e-1 + (0.31750557067975068584e-2 + (0.18369907582308672632e-4 + (0.15761063702089457882e-6 + (0.15577638230480894382e-8 + (0.17663868462699097951e-10 + (0.22126732680711111111e-12 + 0.30273474177737853668e-14 * t) * t) * t) * t) * t) * t) * t;
			}
		case 6: {
				double t = 2*y100 - 13;
				return 0.39317207681134336024e-1 + (0.32504779701937539333e-2 + (0.19354426046513400534e-4 + (0.17081646971321290539e-6 + (0.17485733959327106250e-8 + (0.20593687304921961410e-10 + (0.26917401949155555556e-12 + 0.38562123837725712270e-14 * t) * t) * t) * t) * t) * t) * t;
			}
		case 7: {
				double t = 2*y100 - 15;
				return 0.45896976511367738235e-1 + (0.33300031273110976165e-2 + (0.20423005398039037313e-4 + (0.18567412470376467303e-6 + (0.19718038363586588213e-8 + (0.24175006536781219807e-10 + (0.33059982791466666666e-12 + 0.49756574284439426165e-14 * t) * t) * t) * t) * t) * t) * t;
			}
		case 8: {
				double t = 2*y100 - 17;
				return 0.52640192524848962855e-1 + (0.34139883358846720806e-2 + (0.21586390240603337337e-4 + (0.20247136501568904646e-6 + (0.22348696948197102935e-8 + (0.28597516301950162548e-10 + (0.41045502119111111110e-12 + 0.65151614515238361946e-14 * t) * t) * t) * t) * t) * t) * t;
			}
		case 9: {
				double t = 2*y100 - 19;
				return 0.59556171228656770456e-1 + (0.35028374386648914444e-2 + (0.22857246150998562824e-4 + (0.22156372146525190679e-6 + (0.25474171590893813583e-8 + (0.34122390890697400584e-10 + (0.51593189879111111110e-12 + 0.86775076853908006938e-14 * t) * t) * t) * t) * t) * t) * t;
			}
		case 10: {
				 double t = 2*y100 - 21;
				 return 0.66655089485108212551e-1 + (0.35970095381271285568e-2 + (0.24250626164318672928e-4 + (0.24339561521785040536e-6 + (0.29221990406518411415e-8 + (0.41117013527967776467e-10 + (0.65786450716444444445e-12 + 0.11791885745450623331e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 11: {
				 double t = 2*y100 - 23;
				 return 0.73948106345519174661e-1 + (0.36970297216569341748e-2 + (0.25784588137312868792e-4 + (0.26853012002366752770e-6 + (0.33763958861206729592e-8 + (0.50111549981376976397e-10 + (0.85313857496888888890e-12 + 0.16417079927706899860e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 12: {
				 double t = 2*y100 - 25;
				 return 0.81447508065002963203e-1 + (0.38035026606492705117e-2 + (0.27481027572231851896e-4 + (0.29769200731832331364e-6 + (0.39336816287457655076e-8 + (0.61895471132038157624e-10 + (0.11292303213511111111e-11 + 0.23558532213703884304e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 13: {
				 double t = 2*y100 - 27;
				 return 0.89166884027582716628e-1 + (0.39171301322438946014e-2 + (0.29366827260422311668e-4 + (0.33183204390350724895e-6 + (0.46276006281647330524e-8 + (0.77692631378169813324e-10 + (0.15335153258844444444e-11 + 0.35183103415916026911e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 14: {
				 double t = 2*y100 - 29;
				 return 0.97121342888032322019e-1 + (0.40387340353207909514e-2 + (0.31475490395950776930e-4 + (0.37222714227125135042e-6 + (0.55074373178613809996e-8 + (0.99509175283990337944e-10 + (0.21552645758222222222e-11 + 0.55728651431872687605e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 15: {
				 double t = 2*y100 - 31;
				 return 0.10532778218603311137e0 + (0.41692873614065380607e-2 + (0.33849549774889456984e-4 + (0.42064596193692630143e-6 + (0.66494579697622432987e-8 + (0.13094103581931802337e-9 + (0.31896187409777777778e-11 + 0.97271974184476560742e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 16: {
				 double t = 2*y100 - 33;
				 return 0.11380523107427108222e0 + (0.43099572287871821013e-2 + (0.36544324341565929930e-4 + (0.47965044028581857764e-6 + (0.81819034238463698796e-8 + (0.17934133239549647357e-9 + (0.50956666166186293627e-11 + (0.18850487318190638010e-12 + 0.79697813173519853340e-14 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 17: {
				 double t = 2*y100 - 35;
				 return 0.12257529703447467345e0 + (0.44621675710026986366e-2 + (0.39634304721292440285e-4 + (0.55321553769873381819e-6 + (0.10343619428848520870e-7 + (0.26033830170470368088e-9 + (0.87743837749108025357e-11 + (0.34427092430230063401e-12 + 0.10205506615709843189e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 18: {
				 double t = 2*y100 - 37;
				 return 0.13166276955656699478e0 + (0.46276970481783001803e-2 + (0.43225026380496399310e-4 + (0.64799164020016902656e-6 + (0.13580082794704641782e-7 + (0.39839800853954313927e-9 + (0.14431142411840000000e-10 + 0.42193457308830027541e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 19: {
				 double t = 2*y100 - 39;
				 return 0.14109647869803356475e0 + (0.48088424418545347758e-2 + (0.47474504753352150205e-4 + (0.77509866468724360352e-6 + (0.18536851570794291724e-7 + (0.60146623257887570439e-9 + (0.18533978397305276318e-10 + (0.41033845938901048380e-13 - 0.46160680279304825485e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 20: {
				 double t = 2*y100 - 41;
				 return 0.15091057940548936603e0 + (0.50086864672004685703e-2 + (0.52622482832192230762e-4 + (0.95034664722040355212e-6 + (0.25614261331144718769e-7 + (0.80183196716888606252e-9 + (0.12282524750534352272e-10 + (-0.10531774117332273617e-11 - 0.86157181395039646412e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 21: {
				 double t = 2*y100 - 43;
				 return 0.16114648116017010770e0 + (0.52314661581655369795e-2 + (0.59005534545908331315e-4 + (0.11885518333915387760e-5 + (0.33975801443239949256e-7 + (0.82111547144080388610e-9 + (-0.12357674017312854138e-10 + (-0.24355112256914479176e-11 - 0.75155506863572930844e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 22: {
				 double t = 2*y100 - 45;
				 return 0.17185551279680451144e0 + (0.54829002967599420860e-2 + (0.67013226658738082118e-4 + (0.14897400671425088807e-5 + (0.40690283917126153701e-7 + (0.44060872913473778318e-9 + (-0.52641873433280000000e-10 - 0.30940587864543343124e-11 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 23: {
				 double t = 2*y100 - 47;
				 return 0.18310194559815257381e0 + (0.57701559375966953174e-2 + (0.76948789401735193483e-4 + (0.18227569842290822512e-5 + (0.41092208344387212276e-7 + (-0.44009499965694442143e-9 + (-0.92195414685628803451e-10 + (-0.22657389705721753299e-11 + 0.10004784908106839254e-12 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 24: {
				 double t = 2*y100 - 49;
				 return 0.19496527191546630345e0 + (0.61010853144364724856e-2 + (0.88812881056342004864e-4 + (0.21180686746360261031e-5 + (0.30652145555130049203e-7 + (-0.16841328574105890409e-8 + (-0.11008129460612823934e-9 + (-0.12180794204544515779e-12 + 0.15703325634590334097e-12 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 25: {
				 double t = 2*y100 - 51;
				 return 0.20754006813966575720e0 + (0.64825787724922073908e-2 + (0.10209599627522311893e-3 + (0.22785233392557600468e-5 + (0.73495224449907568402e-8 + (-0.29442705974150112783e-8 + (-0.94082603434315016546e-10 + (0.23609990400179321267e-11 + 0.14141908654269023788e-12 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 26: {
				 double t = 2*y100 - 53;
				 return 0.22093185554845172146e0 + (0.69182878150187964499e-2 + (0.11568723331156335712e-3 + (0.22060577946323627739e-5 + (-0.26929730679360840096e-7 + (-0.38176506152362058013e-8 + (-0.47399503861054459243e-10 + (0.40953700187172127264e-11 + 0.69157730376118511127e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 27: {
				 double t = 2*y100 - 55;
				 return 0.23524827304057813918e0 + (0.74063350762008734520e-2 + (0.12796333874615790348e-3 + (0.18327267316171054273e-5 + (-0.66742910737957100098e-7 + (-0.40204740975496797870e-8 + (0.14515984139495745330e-10 + (0.44921608954536047975e-11 - 0.18583341338983776219e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 28: {
				 double t = 2*y100 - 57;
				 return 0.25058626331812744775e0 + (0.79377285151602061328e-2 + (0.13704268650417478346e-3 + (0.11427511739544695861e-5 + (-0.10485442447768377485e-6 + (-0.34850364756499369763e-8 + (0.72656453829502179208e-10 + (0.36195460197779299406e-11 - 0.84882136022200714710e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 29: {
				 double t = 2*y100 - 59;
				 return 0.26701724900280689785e0 + (0.84959936119625864274e-2 + (0.14112359443938883232e-3 + (0.17800427288596909634e-6 + (-0.13443492107643109071e-6 + (-0.23512456315677680293e-8 + (0.11245846264695936769e-9 + (0.19850501334649565404e-11 - 0.11284666134635050832e-12 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 30: {
				 double t = 2*y100 - 61;
				 return 0.28457293586253654144e0 + (0.90581563892650431899e-2 + (0.13880520331140646738e-3 + (-0.97262302362522896157e-6 + (-0.15077100040254187366e-6 + (-0.88574317464577116689e-9 + (0.12760311125637474581e-9 + (0.20155151018282695055e-12 - 0.10514169375181734921e-12 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 31: {
				 double t = 2*y100 - 63;
				 return 0.30323425595617385705e0 + (0.95968346790597422934e-2 + (0.12931067776725883939e-3 + (-0.21938741702795543986e-5 + (-0.15202888584907373963e-6 + (0.61788350541116331411e-9 + (0.11957835742791248256e-9 + (-0.12598179834007710908e-11 - 0.75151817129574614194e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 32: {
				 double t = 2*y100 - 65;
				 return 0.32292521181517384379e0 + (0.10082957727001199408e-1 + (0.11257589426154962226e-3 + (-0.33670890319327881129e-5 + (-0.13910529040004008158e-6 + (0.19170714373047512945e-8 + (0.94840222377720494290e-10 + (-0.21650018351795353201e-11 - 0.37875211678024922689e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 33: {
				 double t = 2*y100 - 67;
				 return 0.34351233557911753862e0 + (0.10488575435572745309e-1 + (0.89209444197248726614e-4 + (-0.43893459576483345364e-5 + (-0.11488595830450424419e-6 + (0.28599494117122464806e-8 + (0.61537542799857777779e-10 - 0.24935749227658002212e-11 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 34: {
				 double t = 2*y100 - 69;
				 return 0.36480946642143669093e0 + (0.10789304203431861366e-1 + (0.60357993745283076834e-4 + (-0.51855862174130669389e-5 + (-0.83291664087289801313e-7 + (0.33898011178582671546e-8 + (0.27082948188277716482e-10 + (-0.23603379397408694974e-11 + 0.19328087692252869842e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 35: {
				 double t = 2*y100 - 71;
				 return 0.38658679935694939199e0 + (0.10966119158288804999e-1 + (0.27521612041849561426e-4 + (-0.57132774537670953638e-5 + (-0.48404772799207914899e-7 + (0.35268354132474570493e-8 + (-0.32383477652514618094e-11 + (-0.19334202915190442501e-11 + 0.32333189861286460270e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 36: {
				 double t = 2*y100 - 73;
				 return 0.40858275583808707870e0 + (0.11006378016848466550e-1 + (-0.76396376685213286033e-5 + (-0.59609835484245791439e-5 + (-0.13834610033859313213e-7 + (0.33406952974861448790e-8 + (-0.26474915974296612559e-10 + (-0.13750229270354351983e-11 + 0.36169366979417390637e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 37: {
				 double t = 2*y100 - 75;
				 return 0.43051714914006682977e0 + (0.10904106549500816155e-1 + (-0.43477527256787216909e-4 + (-0.59429739547798343948e-5 + (0.17639200194091885949e-7 + (0.29235991689639918688e-8 + (-0.41718791216277812879e-10 + (-0.81023337739508049606e-12 + 0.33618915934461994428e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 38: {
				 double t = 2*y100 - 77;
				 return 0.45210428135559607406e0 + (0.10659670756384400554e-1 + (-0.78488639913256978087e-4 + (-0.56919860886214735936e-5 + (0.44181850467477733407e-7 + (0.23694306174312688151e-8 + (-0.49492621596685443247e-10 + (-0.31827275712126287222e-12 + 0.27494438742721623654e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 39: {
				 double t = 2*y100 - 79;
				 return 0.47306491195005224077e0 + (0.10279006119745977570e-1 + (-0.11140268171830478306e-3 + (-0.52518035247451432069e-5 + (0.64846898158889479518e-7 + (0.17603624837787337662e-8 + (-0.51129481592926104316e-10 + (0.62674584974141049511e-13 + 0.20055478560829935356e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 40: {
				 double t = 2*y100 - 81;
				 return 0.49313638965719857647e0 + (0.97725799114772017662e-2 + (-0.14122854267291533334e-3 + (-0.46707252568834951907e-5 + (0.79421347979319449524e-7 + (0.11603027184324708643e-8 + (-0.48269605844397175946e-10 + (0.32477251431748571219e-12 + 0.12831052634143527985e-13 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 41: {
				 double t = 2*y100 - 83;
				 return 0.51208057433416004042e0 + (0.91542422354009224951e-2 + (-0.16726530230228647275e-3 + (-0.39964621752527649409e-5 + (0.88232252903213171454e-7 + (0.61343113364949928501e-9 + (-0.42516755603130443051e-10 + (0.47910437172240209262e-12 + 0.66784341874437478953e-14 * t) * t) * t) * t) * t) * t) * t) * t;
			 }
		case 42: {
				 double t = 2*y100 - 85;
				 return 0.52968945458607484524e0 + (0.84400880445116786088e-2 + (-0.18908729783854258774e-3 + (-0.32725905467782951931e-5 + (0.91956190588652090659e-7 + (0.14593989152420122909e-9 + (-0.35239490687644444445e-10 + 0.54613829888448694898e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 43: {
				 double t = 2*y100 - 87;
				 return 0.54578857454330070965e0 + (0.76474155195880295311e-2 + (-0.20651230590808213884e-3 + (-0.25364339140543131706e-5 + (0.91455367999510681979e-7 + (-0.23061359005297528898e-9 + (-0.27512928625244444444e-10 + 0.54895806008493285579e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 44: {
				 double t = 2*y100 - 89;
				 return 0.56023851910298493910e0 + (0.67938321739997196804e-2 + (-0.21956066613331411760e-3 + (-0.18181127670443266395e-5 + (0.87650335075416845987e-7 + (-0.51548062050366615977e-9 + (-0.20068462174044444444e-10 + 0.50912654909758187264e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 45: {
				 double t = 2*y100 - 91;
				 return 0.57293478057455721150e0 + (0.58965321010394044087e-2 + (-0.22841145229276575597e-3 + (-0.11404605562013443659e-5 + (0.81430290992322326296e-7 + (-0.71512447242755357629e-9 + (-0.13372664928000000000e-10 + 0.44461498336689298148e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 46: {
				 double t = 2*y100 - 93;
				 return 0.58380635448407827360e0 + (0.49717469530842831182e-2 + (-0.23336001540009645365e-3 + (-0.51952064448608850822e-6 + (0.73596577815411080511e-7 + (-0.84020916763091566035e-9 + (-0.76700972702222222221e-11 + 0.36914462807972467044e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 47: {
				 double t = 2*y100 - 95;
				 return 0.59281340237769489597e0 + (0.40343592069379730568e-2 + (-0.23477963738658326185e-3 + (0.34615944987790224234e-7 + (0.64832803248395814574e-7 + (-0.90329163587627007971e-9 + (-0.30421940400000000000e-11 + 0.29237386653743536669e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 48: {
				 double t = 2*y100 - 97;
				 return 0.59994428743114271918e0 + (0.30976579788271744329e-2 + (-0.23308875765700082835e-3 + (0.51681681023846925160e-6 + (0.55694594264948268169e-7 + (-0.91719117313243464652e-9 + (0.53982743680000000000e-12 + 0.22050829296187771142e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 49: {
				 double t = 2*y100 - 99;
				 return 0.60521224471819875444e0 + (0.21732138012345456060e-2 + (-0.22872428969625997456e-3 + (0.92588959922653404233e-6 + (0.46612665806531930684e-7 + (-0.89393722514414153351e-9 + (0.31718550353777777778e-11 + 0.15705458816080549117e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 50: {
				 double t = 2*y100 - 101;
				 return 0.60865189969791123620e0 + (0.12708480848877451719e-2 + (-0.22212090111534847166e-3 + (0.12636236031532793467e-5 + (0.37904037100232937574e-7 + (-0.84417089968101223519e-9 + (0.49843180828444444445e-11 + 0.10355439441049048273e-12 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 51: {
				 double t = 2*y100 - 103;
				 return 0.61031580103499200191e0 + (0.39867436055861038223e-3 + (-0.21369573439579869291e-3 + (0.15339402129026183670e-5 + (0.29787479206646594442e-7 + (-0.77687792914228632974e-9 + (0.61192452741333333334e-11 + 0.60216691829459295780e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 52: {
				 double t = 2*y100 - 105;
				 return 0.61027109047879835868e0 + (-0.43680904508059878254e-3 + (-0.20383783788303894442e-3 + (0.17421743090883439959e-5 + (0.22400425572175715576e-7 + (-0.69934719320045128997e-9 + (0.67152759655111111110e-11 + 0.26419960042578359995e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 53: {
				 double t = 2*y100 - 107;
				 return 0.60859639489217430521e0 + (-0.12305921390962936873e-2 + (-0.19290150253894682629e-3 + (0.18944904654478310128e-5 + (0.15815530398618149110e-7 + (-0.61726850580964876070e-9 + 0.68987888999111111110e-11 * t) * t) * t) * t) * t) * t;
			 }
		case 54: {
				 double t = 2*y100 - 109;
				 return 0.60537899426486075181e0 + (-0.19790062241395705751e-2 + (-0.18120271393047062253e-3 + (0.19974264162313241405e-5 + (0.10055795094298172492e-7 + (-0.53491997919318263593e-9 + (0.67794550295111111110e-11 - 0.17059208095741511603e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 55: {
				 double t = 2*y100 - 111;
				 return 0.60071229457904110537e0 + (-0.26795676776166354354e-2 + (-0.16901799553627508781e-3 + (0.20575498324332621581e-5 + (0.51077165074461745053e-8 + (-0.45536079828057221858e-9 + (0.64488005516444444445e-11 - 0.29311677573152766338e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 56: {
				 double t = 2*y100 - 113;
				 return 0.59469361520112714738e0 + (-0.33308208190600993470e-2 + (-0.15658501295912405679e-3 + (0.20812116912895417272e-5 + (0.93227468760614182021e-9 + (-0.38066673740116080415e-9 + (0.59806790359111111110e-11 - 0.36887077278950440597e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 57: {
				 double t = 2*y100 - 115;
				 return 0.58742228631775388268e0 + (-0.39321858196059227251e-2 + (-0.14410441141450122535e-3 + (0.20743790018404020716e-5 + (-0.25261903811221913762e-8 + (-0.31212416519526924318e-9 + (0.54328422462222222221e-11 - 0.40864152484979815972e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 58: {
				 double t = 2*y100 - 117;
				 return 0.57899804200033018447e0 + (-0.44838157005618913447e-2 + (-0.13174245966501437965e-3 + (0.20425306888294362674e-5 + (-0.53330296023875447782e-8 + (-0.25041289435539821014e-9 + (0.48490437205333333334e-11 - 0.42162206939169045177e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 59: {
				 double t = 2*y100 - 119;
				 return 0.56951968796931245974e0 + (-0.49864649488074868952e-2 + (-0.11963416583477567125e-3 + (0.19906021780991036425e-5 + (-0.75580140299436494248e-8 + (-0.19576060961919820491e-9 + (0.42613011928888888890e-11 - 0.41539443304115604377e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 60: {
				 double t = 2*y100 - 121;
				 return 0.55908401930063918964e0 + (-0.54413711036826877753e-2 + (-0.10788661102511914628e-3 + (0.19229663322982839331e-5 + (-0.92714731195118129616e-8 + (-0.14807038677197394186e-9 + (0.36920870298666666666e-11 - 0.39603726688419162617e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 61: {
				 double t = 2*y100 - 123;
				 return 0.54778496152925675315e0 + (-0.58501497933213396670e-2 + (-0.96582314317855227421e-4 + (0.18434405235069270228e-5 + (-0.10541580254317078711e-7 + (-0.10702303407788943498e-9 + (0.31563175582222222222e-11 - 0.36829748079110481422e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 62: {
				 double t = 2*y100 - 125;
				 return 0.53571290831682823999e0 + (-0.62147030670760791791e-2 + (-0.85782497917111760790e-4 + (0.17553116363443470478e-5 + (-0.11432547349815541084e-7 + (-0.72157091369041330520e-10 + (0.26630811607111111111e-11 - 0.33578660425893164084e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 63: {
				 double t = 2*y100 - 127;
				 return 0.52295422962048434978e0 + (-0.65371404367776320720e-2 + (-0.75530164941473343780e-4 + (0.16613725797181276790e-5 + (-0.12003521296598910761e-7 + (-0.42929753689181106171e-10 + (0.22170894940444444444e-11 - 0.30117697501065110505e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 64: {
				 double t = 2*y100 - 129;
				 return 0.50959092577577886140e0 + (-0.68197117603118591766e-2 + (-0.65852936198953623307e-4 + (0.15639654113906716939e-5 + (-0.12308007991056524902e-7 + (-0.18761997536910939570e-10 + (0.18198628922666666667e-11 - 0.26638355362285200932e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 65: {
				 double t = 2*y100 - 131;
				 return 0.49570040481823167970e0 + (-0.70647509397614398066e-2 + (-0.56765617728962588218e-4 + (0.14650274449141448497e-5 + (-0.12393681471984051132e-7 + (0.92904351801168955424e-12 + (0.14706755960177777778e-11 - 0.23272455351266325318e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 66: {
				 double t = 2*y100 - 133;
				 return 0.48135536250935238066e0 + (-0.72746293327402359783e-2 + (-0.48272489495730030780e-4 + (0.13661377309113939689e-5 + (-0.12302464447599382189e-7 + (0.16707760028737074907e-10 + (0.11672928324444444444e-11 - 0.20105801424709924499e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 67: {
				 double t = 2*y100 - 135;
				 return 0.46662374675511439448e0 + (-0.74517177649528487002e-2 + (-0.40369318744279128718e-4 + (0.12685621118898535407e-5 + (-0.12070791463315156250e-7 + (0.29105507892605823871e-10 + (0.90653314645333333334e-12 - 0.17189503312102982646e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 68: {
				 double t = 2*y100 - 137;
				 return 0.45156879030168268778e0 + (-0.75983560650033817497e-2 + (-0.33045110380705139759e-4 + (0.11732956732035040896e-5 + (-0.11729986947158201869e-7 + (0.38611905704166441308e-10 + (0.68468768305777777779e-12 - 0.14549134330396754575e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 69: {
				 double t = 2*y100 - 139;
				 return 0.43624909769330896904e0 + (-0.77168291040309554679e-2 + (-0.26283612321339907756e-4 + (0.10811018836893550820e-5 + (-0.11306707563739851552e-7 + (0.45670446788529607380e-10 + (0.49782492549333333334e-12 - 0.12191983967561779442e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 70: {
				 double t = 2*y100 - 141;
				 return 0.42071877443548481181e0 + (-0.78093484015052730097e-2 + (-0.20064596897224934705e-4 + (0.99254806680671890766e-6 + (-0.10823412088884741451e-7 + (0.50677203326904716247e-10 + (0.34200547594666666666e-12 - 0.10112698698356194618e-13 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 71: {
				 double t = 2*y100 - 143;
				 return 0.40502758809710844280e0 + (-0.78780384460872937555e-2 + (-0.14364940764532853112e-4 + (0.90803709228265217384e-6 + (-0.10298832847014466907e-7 + (0.53981671221969478551e-10 + (0.21342751381333333333e-12 - 0.82975901848387729274e-14 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 72: {
				 double t = 2*y100 - 145;
				 return 0.38922115269731446690e0 + (-0.79249269708242064120e-2 + (-0.91595258799106970453e-5 + (0.82783535102217576495e-6 + (-0.97484311059617744437e-8 + (0.55889029041660225629e-10 + (0.10851981336888888889e-12 - 0.67278553237853459757e-14 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 73: {
				 double t = 2*y100 - 147;
				 return 0.37334112915460307335e0 + (-0.79519385109223148791e-2 + (-0.44219833548840469752e-5 + (0.75209719038240314732e-6 + (-0.91848251458553190451e-8 + (0.56663266668051433844e-10 + (0.23995894257777777778e-13 - 0.53819475285389344313e-14 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 74: {
				 double t = 2*y100 - 149;
				 return 0.35742543583374223085e0 + (-0.79608906571527956177e-2 + (-0.12530071050975781198e-6 + (0.68088605744900552505e-6 + (-0.86181844090844164075e-8 + (0.56530784203816176153e-10 + (-0.43120012248888888890e-13 - 0.42372603392496813810e-14 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 75: {
				 double t = 2*y100 - 151;
				 return 0.34150846431979618536e0 + (-0.79534924968773806029e-2 + (0.37576885610891515813e-5 + (0.61419263633090524326e-6 + (-0.80565865409945960125e-8 + (0.55684175248749269411e-10 + (-0.95486860764444444445e-13 - 0.32712946432984510595e-14 * t) * t) * t) * t) * t) * t) * t;
			 }
		case 76: {
				 double t = 2*y100 - 153;
				 return 0.32562129649136346824e0 + (-0.79313448067948884309e-2 + (0.72539159933545300034e-5 + (0.55195028297415503083e-6 + (-0.75063365335570475258e-8 + (0.54281686749699595941e-10 - 0.13545424295111111111e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 77: {
				 double t = 2*y100 - 155;
				 return 0.30979191977078391864e0 + (-0.78959416264207333695e-2 + (0.10389774377677210794e-4 + (0.49404804463196316464e-6 + (-0.69722488229411164685e-8 + (0.52469254655951393842e-10 - 0.16507860650666666667e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 78: {
				 double t = 2*y100 - 157;
				 return 0.29404543811214459904e0 + (-0.78486728990364155356e-2 + (0.13190885683106990459e-4 + (0.44034158861387909694e-6 + (-0.64578942561562616481e-8 + (0.50354306498006928984e-10 - 0.18614473550222222222e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 79: {
				 double t = 2*y100 - 159;
				 return 0.27840427686253660515e0 + (-0.77908279176252742013e-2 + (0.15681928798708548349e-4 + (0.39066226205099807573e-6 + (-0.59658144820660420814e-8 + (0.48030086420373141763e-10 - 0.20018995173333333333e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 80: {
				 double t = 2*y100 - 161;
				 return 0.26288838011163800908e0 + (-0.77235993576119469018e-2 + (0.17886516796198660969e-4 + (0.34482457073472497720e-6 + (-0.54977066551955420066e-8 + (0.45572749379147269213e-10 - 0.20852924954666666667e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 81: {
				 double t = 2*y100 - 163;
				 return 0.24751539954181029717e0 + (-0.76480877165290370975e-2 + (0.19827114835033977049e-4 + (0.30263228619976332110e-6 + (-0.50545814570120129947e-8 + (0.43043879374212005966e-10 - 0.21228012028444444444e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 82: {
				 double t = 2*y100 - 165;
				 return 0.23230087411688914593e0 + (-0.75653060136384041587e-2 + (0.21524991113020016415e-4 + (0.26388338542539382413e-6 + (-0.46368974069671446622e-8 + (0.40492715758206515307e-10 - 0.21238627815111111111e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 83: {
				 double t = 2*y100 - 167;
				 return 0.21725840021297341931e0 + (-0.74761846305979730439e-2 + (0.23000194404129495243e-4 + (0.22837400135642906796e-6 + (-0.42446743058417541277e-8 + (0.37958104071765923728e-10 - 0.20963978568888888889e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 84: {
				 double t = 2*y100 - 169;
				 return 0.20239979200788191491e0 + (-0.73815761980493466516e-2 + (0.24271552727631854013e-4 + (0.19590154043390012843e-6 + (-0.38775884642456551753e-8 + (0.35470192372162901168e-10 - 0.20470131678222222222e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 85: {
				 double t = 2*y100 - 171;
				 return 0.18773523211558098962e0 + (-0.72822604530339834448e-2 + (0.25356688567841293697e-4 + (0.16626710297744290016e-6 + (-0.35350521468015310830e-8 + (0.33051896213898864306e-10 - 0.19811844544000000000e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 86: {
				 double t = 2*y100 - 173;
				 return 0.17327341258479649442e0 + (-0.71789490089142761950e-2 + (0.26272046822383820476e-4 + (0.13927732375657362345e-6 + (-0.32162794266956859603e-8 + (0.30720156036105652035e-10 - 0.19034196304000000000e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 87: {
				 double t = 2*y100 - 175;
				 return 0.15902166648328672043e0 + (-0.70722899934245504034e-2 + (0.27032932310132226025e-4 + (0.11474573347816568279e-6 + (-0.29203404091754665063e-8 + (0.28487010262547971859e-10 - 0.18174029063111111111e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 88: {
				 double t = 2*y100 - 177;
				 return 0.14498609036610283865e0 + (-0.69628725220045029273e-2 + (0.27653554229160596221e-4 + (0.92493727167393036470e-7 + (-0.26462055548683583849e-8 + (0.26360506250989943739e-10 - 0.17261211260444444444e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 89: {
				 double t = 2*y100 - 179;
				 return 0.13117165798208050667e0 + (-0.68512309830281084723e-2 + (0.28147075431133863774e-4 + (0.72351212437979583441e-7 + (-0.23927816200314358570e-8 + (0.24345469651209833155e-10 - 0.16319736960000000000e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 90: {
				 double t = 2*y100 - 181;
				 return 0.11758232561160626306e0 + (-0.67378491192463392927e-2 + (0.28525664781722907847e-4 + (0.54156999310046790024e-7 + (-0.21589405340123827823e-8 + (0.22444150951727334619e-10 - 0.15368675584000000000e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 91: {
				 double t = 2*y100 - 183;
				 return 0.10422112945361673560e0 + (-0.66231638959845581564e-2 + (0.28800551216363918088e-4 + (0.37758983397952149613e-7 + (-0.19435423557038933431e-8 + (0.20656766125421362458e-10 - 0.14422990012444444444e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 92: {
				 double t = 2*y100 - 185;
				 return 0.91090275493541084785e-1 + (-0.65075691516115160062e-2 + (0.28982078385527224867e-4 + (0.23014165807643012781e-7 + (-0.17454532910249875958e-8 + (0.18981946442680092373e-10 - 0.13494234691555555556e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 93: {
				 double t = 2*y100 - 187;
				 return 0.78191222288771379358e-1 + (-0.63914190297303976434e-2 + (0.29079759021299682675e-4 + (0.97885458059415717014e-8 + (-0.15635596116134296819e-8 + (0.17417110744051331974e-10 - 0.12591151763555555556e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 94: {
				 double t = 2*y100 - 189;
				 return 0.65524757106147402224e-1 + (-0.62750311956082444159e-2 + (0.29102328354323449795e-4 + (-0.20430838882727954582e-8 + (-0.13967781903855367270e-8 + (0.15958771833747057569e-10 - 0.11720175765333333333e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 95: {
				 double t = 2*y100 - 191;
				 return 0.53091065838453612773e-1 + (-0.61586898417077043662e-2 + (0.29057796072960100710e-4 + (-0.12597414620517987536e-7 + (-0.12440642607426861943e-8 + (0.14602787128447932137e-10 - 0.10885859114666666667e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 96: {
				 double t = 2*y100 - 193;
				 return 0.40889797115352738582e-1 + (-0.60426484889413678200e-2 + (0.28953496450191694606e-4 + (-0.21982952021823718400e-7 + (-0.11044169117553026211e-8 + (0.13344562332430552171e-10 - 0.10091231402844444444e-12 * t) * t) * t) * t) * t) * t;
			 }
		case 97: case 98:
		case 99: case 100: { // use Taylor expansion for small x (|x| <= 0.0309...)
					   //  (2/sqrt(pi)) * (x - 2/3 x^3  + 4/15 x^5  - 8/105 x^7 + 16/945 x^9) 
					   double x2 = x*x;
					   return x * (1.1283791670955125739
							   - x2 * (0.75225277806367504925
								   - x2 * (0.30090111122547001970
									   - x2 * (0.085971746064420005629
										   - x2 * 0.016931216931216931217))));
				   }
	}
	/* Since 0 <= y100 < 101, this is only reached if x is NaN,
	   in which case we should return NaN. */
	return NAN;
} // w_im_y100


/******************************************************************************/
/*  Library function im_w_of_z                                                */
/******************************************************************************/

	static double
im_w_of_x(double x)
{

	// Steven G. Johnson, October 2012.

	// Uses methods similar to the erfcx calculation:
	// continued fractions for large |x|,
	// a lookup table of Chebyshev polynomials for smaller |x|,
	// and finally a Taylor expansion for |x|<0.01.

	const double ispi = 0.56418958354775628694807945156; // 1 / sqrt(pi)

	if (x >= 0) {
		if (x > 45) { // continued-fraction expansion is faster
			if (x > 5e7) // 1-term expansion, important to avoid overflow
				return ispi / x;
			/* 5-term expansion (rely on compiler for CSE), simplified from:
			   ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
			return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
		}
		return w_im_y100(100/(1+x), x);
	}
	else { // = -im_w_of_x(-x)
		if (x < -45) { // continued-fraction expansion is faster
			if (x < -5e7) // 1-term expansion, important to avoid overflow
				return ispi / x;
			/* 5-term expansion (rely on compiler for CSE), simplified from:
			   ispi / (x-0.5/(x-1/(x-1.5/(x-2/x))))  */
			return ispi*((x*x) * (x*x-4.5) + 2) / (x * ((x*x) * (x*x-5) + 3.75));
		}
		return -w_im_y100(100/(1-x), -x);
	}
} // im_w_of_z

	double
dawson_integral(double x)
{

	// Compute dawson(x) = sqrt(pi)/2 * exp(-x^2) * erfi(x),
	// Dawson's integral for a real argument.

	return spi2 * im_w_of_x(x);
}

/* specfunc/expint.c
 * 
 * Copyright (C) 2007 Brian Gough
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU 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
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

/*
   Chebyshev expansions: based on SLATEC e1.f, W. Fullerton

   Series for AE11       on the interval -1.00000D-01 to  0.
   with weighted error   1.76E-17
   log weighted error  16.75
   significant figures required  15.70
   decimal places required  17.55


   Series for AE12       on the interval -2.50000D-01 to -1.00000D-01
   with weighted error   5.83E-17
   log weighted error  16.23
   significant figures required  15.76
   decimal places required  16.93


   Series for E11        on the interval -4.00000D+00 to -1.00000D+00
   with weighted error   1.08E-18
   log weighted error  17.97
   significant figures required  19.02
   decimal places required  18.61


   Series for E12        on the interval -1.00000D+00 to  1.00000D+00
   with weighted error   3.15E-18
   log weighted error  17.50
   approx significant figures required  15.8
   decimal places required  18.10


   Series for AE13       on the interval  2.50000D-01 to  1.00000D+00
   with weighted error   2.34E-17
   log weighted error  16.63
   significant figures required  16.14
   decimal places required  17.33


   Series for AE14       on the interval  0.          to  2.50000D-01
   with weighted error   5.41E-17
   log weighted error  16.27
   significant figures required  15.38
   decimal places required  16.97
 */

static const ChebyshevSeries<double, 39> E1a1_series = {
	38, 38, 20, {
		0.121503239716065790,
		-0.065088778513550150,
		0.004897651357459670,
		-0.000649237843027216,
		0.000093840434587471,
		0.000000420236380882,
		-0.000008113374735904,
		0.000002804247688663,
		0.000000056487164441,
		-0.000000344809174450,
		0.000000058209273578,
		0.000000038711426349,
		-0.000000012453235014,
		-0.000000005118504888,
		0.000000002148771527,
		0.000000000868459898,
		-0.000000000343650105,
		-0.000000000179796603,
		0.000000000047442060,
		0.000000000040423282,
		-0.000000000003543928,
		-0.000000000008853444,
		-0.000000000000960151,
		0.000000000001692921,
		0.000000000000607990,
		-0.000000000000224338,
		-0.000000000000200327,
		-0.000000000000006246,
		0.000000000000045571,
		0.000000000000016383,
		-0.000000000000005561,
		-0.000000000000006074,
		-0.000000000000000862,
		0.000000000000001223,
		0.000000000000000716,
		-0.000000000000000024,
		-0.000000000000000201,
		-0.000000000000000082,
		0.000000000000000017,
	}
};

static const ChebyshevSeries<double, 25> E1a2_series = {
	24, 24, 15,
	{
		0.582417495134726740,
		-0.158348850905782750,
		-0.006764275590323141,
		0.005125843950185725,
		0.000435232492169391,
		-0.000143613366305483,
		-0.000041801320556301,
		-0.000002713395758640,
		0.000001151381913647,
		0.000000420650022012,
		0.000000066581901391,
		0.000000000662143777,
		-0.000000002844104870,
		-0.000000000940724197,
		-0.000000000177476602,
		-0.000000000015830222,
		0.000000000002905732,
		0.000000000001769356,
		0.000000000000492735,
		0.000000000000093709,
		0.000000000000010707,
		-0.000000000000000537,
		-0.000000000000000716,
		-0.000000000000000244,
		-0.000000000000000058,
	}
};

static const ChebyshevSeries<double, 19> E11_series = {
	18, 18, 13,
	{
		-16.11346165557149402600,
		7.79407277874268027690,
		-1.95540581886314195070,
		0.37337293866277945612,
		-0.05692503191092901938,
		0.00721107776966009185,
		-0.00078104901449841593,
		0.00007388093356262168,
		-0.00000620286187580820,
		0.00000046816002303176,
		-0.00000003209288853329,
		0.00000000201519974874,
		-0.00000000011673686816,
		0.00000000000627627066,
		-0.00000000000031481541,
		0.00000000000001479904,
		-0.00000000000000065457,
		0.00000000000000002733,
		-0.00000000000000000108
	}
};

static const ChebyshevSeries<double, 16> E12_series = {
	15, 15, 10,
	{
		-0.03739021479220279500,
		0.04272398606220957700,
		-0.13031820798497005440,
		0.01441912402469889073,
		-0.00134617078051068022,
		0.00010731029253063780,
		-0.00000742999951611943,
		0.00000045377325690753,
		-0.00000002476417211390,
		0.00000000122076581374,
		-0.00000000005485141480,
		0.00000000000226362142,
		-0.00000000000008635897,
		0.00000000000000306291,
		-0.00000000000000010148,
		0.00000000000000000315
	}
};

static const ChebyshevSeries<double, 25> E13_series = {
	24, 24, 15,
	{
		-0.605773246640603460,
		-0.112535243483660900,
		0.013432266247902779,
		-0.001926845187381145,
		0.000309118337720603,
		-0.000053564132129618,
		0.000009827812880247,
		-0.000001885368984916,
		0.000000374943193568,
		-0.000000076823455870,
		0.000000016143270567,
		-0.000000003466802211,
		0.000000000758754209,
		-0.000000000168864333,
		0.000000000038145706,
		-0.000000000008733026,
		0.000000000002023672,
		-0.000000000000474132,
		0.000000000000112211,
		-0.000000000000026804,
		0.000000000000006457,
		-0.000000000000001568,
		0.000000000000000383,
		-0.000000000000000094,
		0.000000000000000023
	}
};

static const ChebyshevSeries<double, 26> E14_series = {
	25, 25, 13,
	{
		-0.18929180007530170,
		-0.08648117855259871,
		0.00722410154374659,
		-0.00080975594575573,
		0.00010999134432661,
		-0.00001717332998937,
		0.00000298562751447,
		-0.00000056596491457,
		0.00000011526808397,
		-0.00000002495030440,
		0.00000000569232420,
		-0.00000000135995766,
		0.00000000033846628,
		-0.00000000008737853,
		0.00000000002331588,
		-0.00000000000641148,
		0.00000000000181224,
		-0.00000000000052538,
		0.00000000000015592,
		-0.00000000000004729,
		0.00000000000001463,
		-0.00000000000000461,
		0.00000000000000148,
		-0.00000000000000048,
		0.00000000000000016,
		-0.00000000000000005
	}
};

	static double
E1_scaled(const double x)
{
	if (x <= -10.0)
		return (1.0 + E1a1_series.eval_double(20.0/x + 1.0))/x;
	else if (x <= -4.0)
		return (1.0 + E1a2_series.eval_double((40.0/x + 7.0)/3.0))/x;
	else if (x <= -1.0)
		return exp(x) * (-log(fabs(x))
				+ E11_series.eval_double((2.0*x + 5.0)/3.0));
	else if (x == 0.0)
		return NAN;
	else if (x <= 1.0)
		return exp(x) * (-log(fabs(x)) - 0.6875 + x
				+ E12_series.eval_double(x));
	else if (x <= 4.0)
		return (1.0 + E13_series.eval_double((8.0/x - 5.0)/3.0))/x;
	else
		return (1.0 +  E14_series.eval_double(8.0/x-1.0))/x;
}

	double
Ei_scaled(double x)
{
	return -E1_scaled(-x);
}

	long double
Ei(long double x)
{
	return expl(x)*Ei_scaled(x);
}

Complex epsIBT3(double E,double N, double Eg,double Eh)
{
	long double Eg2=Eg*Eg;
	long double Eh2=Eh*Eh;
	long double E2=E*E;
	long double aEmEg=fabs(E)-Eg;
	long double aEmEh=fabs(E)-Eh;
	long double logEhEg=logl(Eh/Eg);
	long double Eg2Eh2=Eg2*Eh2;
	long double Eh2mEg2=Eh2-Eg2;
	long double C=(Eh2mEg2)*(Eh2+Eg2-8*Eg*Eh)/12+Eg2Eh2*logEhEg;
	long double N_C=N/C;
	long double N_CE2=N_C/E2;

	double ei=( fabs(E)>Eg && fabs(E)<Eh ? (E>0?1:-1)*N_CE2*aEmEg*aEmEg*aEmEh*aEmEh : 0.0 );

	long double x=2*E*(Eh*(Eg2+E2)+Eg*(Eh2+E2));
	long double y=E2*(Eh2+Eg2+4*Eg*Eh+E2)+Eg2Eh2;
	long double a=y-x;
	long double b=y+x;
	long double c=-2*Eg2Eh2*logEhEg;
	long double d=-3*Eh2mEg2;
	long double N_piCE2=N_CE2/M_PI;
	long double N_piC=N_C/M_PI;

	double er;
	if(fabs(E)<2e-4*Eg) er=N_piC*(2*(Eh2+Eg2+4*Eg*Eh)*logEhEg+2*d);
	else if(E==Eg || E==Eh) er=N_piCE2*(b*logl(fabsl((Eh+E)/(Eg+E)))+c+d*E2);
	else if(E==-Eg || E==-Eh) er=N_piCE2*(a*logl(fabsl((Eh-E)/(Eg-E)))+c+d*E2);
	else er=N_piCE2*(a*logl(fabsl((Eh-E)/(Eg-E)))+b*logl(fabsl((Eh+E)/(Eg+E)))+c+d*E2);

	return Complex(er,ei);
}

Complex epsIBTL5(double E,double N, double Eg,double Eh,double Ec,double B)
{
	long double Eg2=Eg*Eg;
	long double Eh2=Eh*Eh;
	long double Ec2=Ec*Ec;
	long double E2=E*E;
	long double B2=B*B;
	long double aEmEg=fabs(E)-Eg;
	long double aEmEh=fabs(E)-Eh;
	long double aEmEc=fabs(E)-Ec;
	long double logEhEg=logl(Eh/Eg);
	long double EgEh=Eg*Eh;
	long double Eg2Eh2=Eg2*Eh2;
	long double EhmEg=Eh-Eg;
	long double alpha=Ec2+B2;
	long double alpha2=alpha*alpha;
	long double beta=Eg+Eh;
	long double gamma=Eg2+Eh2;
	long double delta=alpha2*(beta-Ec)-alpha*beta*EgEh+Eg2Eh2*Ec;
	long double eta=alpha2-alpha*gamma+4*alpha*(beta*Ec-Ec2-EgEh)+Eg2Eh2;
	long double EcmEg=Ec-Eg;
	long double EhmEc=Eh-Ec;
	long double A=atanl((EcmEg)/B)+atanl((EhmEc)/B);
	long double L=logl((EhmEc*EhmEc+B2)/(EcmEg*EcmEg+B2));
	long double C=(EhmEg)*(4*Ec-3*beta)/2+Eg2Eh2/alpha*logEhEg+(4*delta*A-eta*(B*L+2*Ec*A))/(2*alpha*B);
	long double N_C=N/C;
	long double N_CE2=N_C/E2;

	double ei=( fabs(E)>Eg && fabs(E)<Eh ? (E>0?1:-1)*N_CE2*aEmEg*aEmEg*aEmEh*aEmEh/(aEmEc*aEmEc+B2) : 0.0 );

	long double x=2*beta*E*(E2+EgEh);
	long double y=E2*(E2+gamma+4*EgEh)+Eg2Eh2;
	long double a=(y-x)/(E2-2*Ec*E+alpha);
	long double b=(y+x)/(E2+2*Ec*E+alpha);
	long double d=delta*E2+alpha2*(alpha*beta-beta*EgEh-gamma*Ec)+4*alpha*EgEh*Ec*(beta*Ec-alpha)+Eg2Eh2*(2*alpha*Ec-4*Ec*Ec2);
	long double e=eta*E2+alpha2*(alpha-gamma)-4*alpha2*EgEh+alpha*Eg2Eh2+4*alpha*beta*Ec*EgEh-4*Ec2*Eg2Eh2;;
	long double E2palpha=E2+alpha;
	long double c0=(4*A*d-(B*L+2*Ec*A)*e)/(alpha*B*(4*Ec2*E2-E2palpha*E2palpha));
	long double c=c0*E2-2*Eg2Eh2/alpha*logEhEg;
	long double N_piC=N_C/M_PI;
	long double N_piCE2=N_piC/E2;

	double er;
	if(fabs(E)<2e-5*Eg) er=N_piC*((Eh2-Eg2-4*beta*EhmEg)/alpha+EhmEg*4*EgEh*Ec/alpha2+((8*EgEh+2*gamma)/alpha-(2*Eg2Eh2+8*beta*EgEh*Ec)/alpha2+8*Eg2Eh2*Ec2/alpha/alpha2)*logEhEg+c0);
	else if(E==Eg || E==Eh) er=N_piCE2*(b*logl(fabsl((Eh+E)/(Eg+E)))+c);
	else if(E==-Eg || E==-Eh) er=N_piCE2*(a*logl(fabsl((Eh-E)/(Eg-E)))+c);
	else er=N_piCE2*(a*logl(fabsl((Eh-E)/(Eg-E)))+b*logl(fabsl((Eh+E)/(Eg+E)))+c);

	return Complex(er,ei);
}

Complex epsGP3(double E, double N, double Ec, double B)
{
	if (!N) return 0.0;

	double S = B*M_SQRT2;

	double um = (E - Ec)/S;
	double up = (E + Ec)/S;
	double g = (exp(-um*um) - exp(-up*up))/M_SQRTPI;
	double d = ( dawson_integral(up) - dawson_integral(um))*M_2_PI;
	return Complex(d, g)*N/(S*Ec);
}

Complex epsHET2(double E,double N, double Eg)
{
	double aEmEg=fabs(E)-Eg;
	double E2=E*E;
	double E3=E*E2;
	double N3EgE2=N*3*Eg/E2;
	double ei=( fabs(E)>Eg ? N3EgE2*aEmEg*aEmEg/E3 : 0.0 );

	long double EmEg=E-Eg;
	long double EpEg=E+Eg;
	long double a=-EmEg*EmEg/E3;
	long double b=EpEg*EpEg/E3;
	long double c=-2/(3*Eg);
	long double d=-2*Eg;

	double er;
	if( fabs(E)<Eg*5e-3 ) er=N/(5*M_PI*Eg*Eg);
	else if( E == Eg ) er=N3EgE2/M_PI*(b*logl(2.0)+c+d/E2);
	else if( E == -Eg ) er=N3EgE2/M_PI*(a*logl(2.0)+c+d/E2);
	else if( fabs(E)<Eg ) er=N3EgE2/M_PI*((a*log1pl(-E/Eg)+b*log1pl(E/Eg))+c+d/E2);
	else er=N3EgE2/M_PI*((a*logl(fabsl(1-E/Eg))+b*logl(fabsl(1+E/Eg)))+c+d/E2);

	return Complex(er,ei);
}

Complex epsUT4(double E,double N,double Eu,double Eg,double Eh)
{
	long double Eg_2=Eg/2;
	long double EhpEg_2=Eh+Eg_2;
	long double E2=E*E;
	long double Eg_Eu=Eg/Eu;
	long double Eg_2Eu=Eg_Eu/2;
	long double Eh_Eu=Eh/Eu;
	long double Eh_2Eu=Eh_Eu/2;
	long double EgEh=Eg*Eh;
	long double C=expl(-Eg_Eu)*(Ei(Eg_Eu)-Ei(Eg_2Eu))
		+expl(Eh_Eu)*(Ei(-EhpEg_2/Eu)-Ei(-Eh_Eu))
		+expl(-Eg_2Eu)*logl(Eg_2/EhpEg_2)
		+logl(Eh/Eg)*(1-EgEh/(Eu*(Eh-Eg)))+Eh_2Eu+Eg_2Eu;
	long double NCE2=N/(C*E2);
	long double EhmaEEu=(Eh-fabs(E))/Eu;
	long double aEmEg_Eu=(fabs(E)-Eg)/Eu;
	long double EhmEg=Eh-Eg;
	long double EhpEg=Eh+Eg;
	long double mEhmEgp2aE=2*fabs(E)-Eh-Eg;
	long double Eu2=Eu*Eu;
	long double Eg2=Eg*Eg;
	long double Eh2=Eh*Eh;
	long double EgEu=Eg*Eu;
	long double EhEu=Eh*Eu;

	double ei;
	if(fabs(E)>Eg_2 && fabs(E)<Eg) ei=NCE2*(expl(aEmEg_Eu)-expl(-Eg_2Eu));
	else if(fabs(E)>=Eg && fabs(E)<=Eh) ei=NCE2*(1+(EhmEg)/(4*Eu)-mEhmEgp2aE*mEhmEgp2aE/(4*Eu*EhmEg)-expl(-Eg_2Eu));
	else if(fabs(E)>Eh && fabs(E)<EhpEg_2) ei=NCE2*(expl(EhmaEEu)-expl(-Eg_2Eu));
	else ei=0;
	ei=(E>=0?1:-1)*ei;

	long double NpiCE2=NCE2/M_PI;
	long double x=-E*(EhpEg)/(Eu*(EhmEg));
	long double y=1-(EgEh+E2)/(Eu*(EhmEg));
	long double a=y-x;
	long double b=y+x;
	long double c=2*((EgEh/(Eu*(EhmEg))-1)*logl(Eh/Eg)-Ei_scaled(Eg_Eu)+Ei_scaled(-Eh_Eu)
			+expl(-Eg_2Eu)*(Ei_scaled(Eg_2Eu)-Ei_scaled(-EhpEg_2/Eu)+logl(2*Eh/Eg+1)));

	double er=0.0;
	if(fabs(E)<1e-3) er=N/(M_PI*C)*(EhpEg/(Eh*Eg*Eu)+(Eh*Eh-Eg*Eg)/(Eg2*Eh2)-2/(Eu*EhmEg)*logl(Eh/Eg)
			+1/Eu2*Ei_scaled(Eg_Eu)-1/EgEu-1/Eg2
			-1/Eu2*Ei_scaled(-Eh_Eu)-1/EhEu+1/Eh2
			-expl(-Eg_2Eu)*(16*Eh*EhpEg/(Eg*Eg*(4*Eh*Eh+4*Eh*Eg+Eg*Eg))
				+1/Eu2*Ei_scaled(Eg_2Eu)-2/EgEu-4/Eg2
				-1/Eu2*Ei_scaled(-EhpEg_2/Eu)-1/(EhpEg_2*Eu)+1/(EhpEg_2*EhpEg_2)));
	else {
		if(E==Eg) er+=logl((Eh-Eg)/Eu)+M_EULER-Ei_scaled(-(Eh-E)/Eu);
		else if(E==Eh) er+=logl(Eu/(Eh-Eg))+M_EULER+Ei_scaled((Eg-E)/Eu);
		else er+=a*logl(fabs((Eh-E)/(Eg-E)))+Ei_scaled((Eg-E)/Eu)-Ei_scaled(-(Eh-E)/Eu);
		if(E==-Eg) er+=logl((Eh-Eg)/Eu)+M_EULER-Ei_scaled(-(Eh+E)/Eu);
		else if(E==-Eh) er+=logl(Eu/(Eh-Eg))+M_EULER+Ei_scaled((Eg+E)/Eu);
		else er+=b*logl(fabs((Eh+E)/(Eg+E)))+Ei_scaled((Eg+E)/Eu)-Ei_scaled(-(Eh+E)/Eu);
		er+=c;
		if(E==Eg_2) er+=-expl(-Eg_2Eu)*(logl(Eh_Eu)+M_EULER-Ei_scaled((E-EhpEg_2)/Eu));
		else if(E==EhpEg_2) er+=-expl(-Eg_2Eu)*(-logl(Eh_Eu)-M_EULER+Ei_scaled((Eg_2-E)/Eu));
		else er+=-expl(-Eg_2Eu)*(logl(fabs(Eh/(Eg_2-E)+1))+Ei_scaled((Eg_2-E)/Eu)-Ei_scaled((E-EhpEg_2)/Eu));
		if(E==-Eg_2) er+=-expl(-Eg_2Eu)*(logl(Eh_Eu)+M_EULER-Ei_scaled(-(E+EhpEg_2)/Eu));
		else if(E==-EhpEg_2)er+=-expl(-Eg_2Eu)*(-logl(Eh_Eu)-M_EULER+Ei_scaled((Eg_2+E)/Eu));
		else er+=-expl(-Eg_2Eu)*(logl(fabs(Eh/(Eg_2+E)+1))+Ei_scaled((Eg_2+E)/Eu)-Ei_scaled(-(E+EhpEg_2)/Eu));
		er=NpiCE2*er;
	}
	return Complex(er,ei);
}

const double kB=8.61734318378e-05;              // = kB/e = 1.3806505e-23/1.60217653e-19

static double fBE(double E,double T)
{
        if(E==0.0) return 0.0;
        double t=fabs(E)/(kB*T);
        if(t>200) return 0.0;
        return 1/(exp(t)-1);
}

Complex udm::df(double E)
{
	const double FWHM_to_rms=1/(2*sqrt(2*M_LN2));
	Complex eps=1;

	for(unsigned j=0;j<nph;j++) eps+=epsGP3(E,Nph[j],tnu_to_E*nuph[j],tnu_to_E*FWHM_to_rms*betaph[j]);

	if(he) eps+=epsHET2(E,Nvx,Ex+Eg/2);

	if(ut) eps+=epsUT4(E,Nut,Eu,Eg,Eh);

	for(unsigned j=0;j<nloc;j++) eps+=epsGP3(E,Nloc[j],Eloc[j],Bloc[j]);

	double sumA=A0;
	for(unsigned j=0;j<nex;j++) sumA+=Aex[j];
	double np=0.0;
	if(dEg) {
		if(Ep) np=fBE(Ep,300);
	}
	eps+=epsIBT3(E,(np+1)/(2*np+1)*A0*Nvc/sumA,Eg+Ep,Eh);
	if(dEg) eps+=epsIBT3(E,np/(2*np+1)*A0*Nvc/sumA,Eg-Ep,Eh);
	if(gex) {
		for(unsigned j=0;j<nex;j++) {
			eps+=epsIBTL5(E,(np+1)/(2*np+1)*(1-Gex[j])*Aex[j]*Nvc/sumA,Eg+Ep,Eh,Eex[j],Bex[j]);
			if(dEg) eps+=epsIBTL5(E,np/(2*np+1)*(1-Gex[j])*Aex[j]*Nvc/sumA,Eg-Ep,Eh,Eex[j],Bex[j]);
			eps+=epsGP3(E,Gex[j]*Aex[j]*Nvc/sumA,Eex[j],Bex[j]/M_SQRT2);
		}
	}
	else {
		for(unsigned j=0;j<nex;j++) {
			eps+=epsIBTL5(E,(np+1)/(2*np+1)*Aex[j]*Nvc/sumA,Eg+Ep,Eh,Eex[j],Bex[j]);
			if(dEg) eps+=epsIBTL5(E,np/(2*np+1)*Aex[j]*Nvc/sumA,Eg-Ep,Eh,Eex[j],Bex[j]);
		}
	}

	return eps;
}

int main(int argc, char *argv[])
{
	try {
		terminate_program err;
		if(argc==1) throw err;
		if(argc!=2 && argc!=3) {
			err+="wrong number of arguments";
			throw err;
		}
		
		string parfile=argv[1];
	
		if(parfile == "--arrange") {
			if(argc==2) throw err;
			
			parfile=argv[2];
			
			udm model(parfile);
			
			model.arrange();
			model.print();
			
			exit(0);
		}

		udm model(parfile);

		double val1=0.01;
		double val2=100.0;
		unsigned npoints=4001;
		string vars="E";
		string dist="exp";

		if(argc==3) {
			string defvar(argv[2]);
			size_t p=defvar.find('=');
			if(p==string::npos) {
				err+="wrong definition of variable";
				throw err;
			}

			vars=defvar.substr(0,p);
			trim(vars);
			if(vars!="E" && vars!="lambda" && vars!="tnu" && vars!="nu") {
				err+="unknown variable "+vars;
				throw err;
			}

			string val1s=defvar.substr(p+1);
			string val2s;
			string npointss;

			p=val1s.find('~');
			if(p!=string::npos) {
				val2s=val1s.substr(p+1);
				trim(val2s);

				val1s=val1s.substr(0,p);
				trim(val1s);

				p=val2s.find(':');
				if(p==string::npos) {
					err+="wrong definition of variable";
					throw err;
				}
				npointss=val2s.substr(p+1);
				trim(npointss);
				val2s=val2s.substr(0,p);
				trim(val2s);

				p=npointss.find(':');
				if(p==string::npos) dist="lin";
				else {
					dist=npointss.substr(p+1);
					trim(dist);
					npointss=npointss.substr(0,p);
					trim(npointss);
				}
			}

			try{ val1=lexical_cast<double>(val1s); }
			catch(bad_lexical_cast) {
				err+="can't covert \""+val1s+"\" to number";
				throw err;
			}
			try{ val2=lexical_cast<double>(val2s); }
			catch(bad_lexical_cast) {
				err+="can't covert \""+val2s+"\" to number";
				throw err;
			}
			try{ npoints=lexical_cast<unsigned>(npointss); }
			catch(bad_lexical_cast) {
				err+="can't covert \""+npointss+"\" to number";
				throw err;
			}
			if(val1<0 || val2<=val1) {
				err+="values of "+vars+" must be positive and value2>value1";
				throw err;
			}
			if(npoints<2) {
				err+="wrong number of points (number_of_points>1)";
				throw err;
			}
			if(dist!="lin" && dist!="exp") {
				err+="wrong distribution ('lin' or 'exp')";
				throw err;
			}
		}

		for(unsigned i=0;i<npoints;i++) {
			double val=val1;
			if(val==0.0 && vars=="lambda") printf("%e\t%e\t%e\n",0.0,1.0,0.0);
			else {
				if(dist=="lin") val+=i*(val2-val1)/(npoints-1);
				if(dist=="exp") val*=exp(i*log(val2/val1)/(npoints-1));
				double E=val;
				if(vars=="lambda") E=lambda_E/val;
				else if(vars=="tnu") E=tnu_to_E*val;
				else if(vars=="nu") E=nu_to_E*val;

				Complex oc=sqrt(model.df(E));
				printf("%e\t%e\t%e\n",val,real(oc),imag(oc));
			}
		}

	} catch(terminate_program err) {
		if( err.ifmsg() ) cerr << err.read() << endl;
		cerr << "Program for calculation of optical constants using Universal Dispersion Model" << endl;
		cerr << "for format of <par-file> see newAD.physics.muni.cz\n" << endl;
		cerr << "author: Daniel Franta, 2015\n" << endl;
		cerr << "usage:\n" << argv[0] << " <par-file> [quantity=value1[~value2:number_of_points[:distribution]]]\n" << endl;
		cerr << "quantity:" << endl;
		cerr << "   E ... photon energy (eV)" << endl;
		cerr << "   lambda ... wavelength (nm)" << endl;
		cerr << "   tnu ... wavenumber (1/cm)" << endl;
		cerr << "   nu ... frequency (THz)\n" << endl;
		cerr << "value1 >= 0" << endl;
		cerr << "value2 > value1" << endl;
		cerr << "number_of_points > 1\n" << endl;
		cerr << "distribution:" << endl;
		cerr << "   lin ... linear distribution of values (default)" << endl;
		cerr << "   exp ... linear distribution of values\n" << endl;
		cerr << "default: E=0.01~100:4001:exp\n" << endl;
		cerr << "usage:\n" << argv[0] << " --arrange <par-file>\n" << endl;
		cerr << "arrange numbered parameters in <par-file> to stdout\n" << endl;
		
		exit(1);
	}
}