/**
* Computer for amplitude reflection coefficients of an electromagnetic
* wave impinging on a coated substrate, where coating and substrate dielectric
* functions (functions of wavelength) are known.
*
* The model is application of Snell's law to compute angles relative to
* the surface normal (the angle in front of the coating to be provided by
* the user). There is a "coherent" and an "incoherent" option of the model
* (only the first relevant to the computation of beam splitters): in the coherent
* computation, the wave is reflected off the air-coating interface, and the
* portion that is transmitted into the coating is travelling twice through
* it (with Fresnel amplitude reflection/transmission coefficients applied at
* each interface, and the phase factor of the passage following from the
* OPD within), 4 times through it etc, all adding up to another component
* after one additional re-transmission to the front surface, and mingling
* with the direct reflection. The product of all these complex Fresnel
* coefficients and phase factors (including a phase factor from the free
* path length through air of the directly reflected wave) defines an amplitude
* reflection coefficient. All these computations are done separately for
* the p and s polarization. The ratio of the two complex amplitude reflection
* coefficients defines another complex number, the argument of which is
* the relative phase shift between the two polarizations.
*
* The similar computation is done for the wave that is transmitted through
* the coating into the substrate, with the difference that no interference
* calculation is done in the substrate. This means a wave can run through
* all three interface and exit at the rear surface, or it can be reflected
* within the coating once or more often and enter the substrate and exit
* through the rear surface, but there is no "coherent" superposition of
* multi-travels through the substrate. As a side effect, the sum of the
* transmission and reflection energies does not add to 1, because the
* light that is entering the substrate may be reflected from the rear
* interface to the air with some reflection coefficient; this is a portion
* that remains trapped in the substrate in the current calculation and
* does not show in transmission or reflection. This means, the effect
* of an anti-reflection coating at the rear surface is not taken into account.
*
* The user command line interface is described near the main() function
* below.
*
* Compilation:
*   make MatisseBS
* 
* Richard J. Mathar www.strw.leidenuniv.nl/~mathar
* 2008-05-29
*
* Literature:
*
* The wavelength dependent reflection of the MIDI N-band beam splitter
* on the coated side is plotted in Fig 3 in the article of Bakker, Przygodda,
* Chesneau, Jaffe in ESA SP-522 "VLTI MIDI scientific observation procedures:
* Synergie between MIDI and GENIE - the MIDI nuller", Proc. GENIE-DARWIN
* workshop - Hungint for Planets, 3-6 June 2002, Leiden.
*
* Phase shifters are discussed by Rabbbia, Gay, Rivet and Schneider in
* "Achromatic Phase shifters: The "mirror" approaches", same proceedings
* as above.
*
* Serabyn and Colavita, Fully symmetric nulling beam combinsers
*   Appl. Opt vol 40  no 10 (2001) 1668
**/
#include <complex>
#include <vector>
#include <algorithm>
#include <iostream>
#include <unistd.h>

#include "asin.cpp"
#include "units.h"

using namespace std ;

/* default substrate thickness: 1 mm */
#define MATISSE_SUBSTR_THICK 1.e-3

/** A complex value along (as a function of) some 1-dim axis.
* @author Richard J. Mathar
* @since 2008-05-27
*/
class Clist {
public:
	/* The coordinate on the axis.
	*/
	double x ;

	/* The complex value at that coordinate
	*/
	complex<double> val ;

	/* Constructor with the pair of values
	*/
	Clist(const double xval, const complex<double> yval) : x(xval), val(yval)
	{
	}
protected:
private:
} ;

/** Report a value to the output stream
* @param[in,out] os the output stream
* @param[in] v the value to be printed
* @return the repositioned output stream
* @author Richard J. Mathar
* @since 2008-05-27
*/
ostream & operator<< (ostream &os, const Clist & v)
{
	os << v.val << " at " << v.x ;
	return os ;
}

/** Comparison of two dielectric values by location on the wavenumber axis
* @param[in] left the dielectric function at the left of the comparison
* @param[in] right the dielectric function at the right of the comparison
* @return true if the left wavenumber is smaller than the right wavenumber
* @author Richard J. Mathar
* @since 2008-05-27
*/
bool operator< (const Clist &left, const Clist & right)
{
	return ( left.x < right.x ) ;
}

/** Dielectric Function.
* @author Richard J. Mathar
* @since 2008-05-27
*/
class DielF {
public:

	/* A table of wavenumbers and associated complex valued
	* dielectric functions at these points.
	*/
	vector<Clist> tab ;

	/* Default Constructor.
	* The list of values is empty by default.
	*/
	DielF()
	{
	}

	/* Constructor
	* @param[in] dfun a table of dielectric functions
	*/
	DielF(const vector<Clist> & dfun) : tab(dfun)
	{
		sort(tab.begin(),tab.end()) ;
	}

	/* Add one point to the dispersion curve
	* @param[in] val the new wavenumber and complex dielectric value
	*/
	void disp(const Clist & val)
	{
		tab.push_back(val) ;
		sort(tab.begin(),tab.end()) ;
	}

	/* Return the value at a particular wavenumber
	* @param[in] wn the vacuum wavenumber [1/cm]
	* @return the value of the bulk dielectric function
	* @warning this uses a simple polynomial interpolation between
	*   all values and ought not be used with more than 4 abscissa points.
	*/
	complex<double> n(const double wn) const
	{
		/* Simple Lagrange interpolation between all values
		*/
		complex<double> rindx(0.) ;
		for(int i=0; i< tab.size() ;i++)
		{
			double p=1. ;
			for(int j=0; j < tab.size() ; j++)
				if ( j != i)
					p *= (wn - tab[j].x)/(tab[i].x-tab[j].x) ;
			rindx += p*tab[i].val ;
		}
		return rindx ;
	}

	/** Fresnel amplitude reflection coefficient.
	* @param[in] phi angle of incidence in this materials side [rad]
	* @param[in] wn the vacuum wavenumber of the elm wave
	* @param[in] oth the other material from which we reflect
	* @param[out] R the p polarixation (index 0) and s polarization (index 1) Fresnel coeffs
	*/
	void fresn_coefR(const double phi, const double wn, const DielF & oth, complex<double> R[2]) const
	{
		/* the relative dielectric constant is the ratio of
		* the other materials dielectric constant over the one here.
		*/
		const complex<double> nratio = oth.n(wn)/n(wn) ;

		/* Handle the case of perpendicular incidence separately to avoid
		* singularities in the formalism.
		*/
		if (phi == 0.)
		{
			R[0] = (nratio-1.)/(nratio+1.) ;
			R[1] = (1.-nratio)/(nratio+1.) ;
		}
		else
		{
			/* Snell's law in the complex version.
			*/
			const complex<double> psi = asin( sin(phi)/nratio) ;
			R[0] = tan(phi-psi)/tan(phi+psi) ;
			R[1] = -sin(phi-psi)/sin(phi+psi) ;
		}
	}

	/** Angle of transmission (relative to surface normal)
	* @param[in] phi angle of incidence [rad]
	* @param[in] wn vacuum wavenumber
	* @param[in] oth the dielectric to transmit into
	* @return the angle inside the other material [rad]
	*/
	double Tang(const double phi, const double wn, const DielF & oth) const
	{
		/* the relative dielectric constant is the ratio of
		* the other materials dielectric constant divided through the one here.
		*/
		const complex<double> nratio = oth.n(wn)/n(wn) ;

		/* Snell's law in the complex version.
		*/
		const complex<double> psi = asin( sin(phi)/nratio) ;
		return psi.real() ;
	}

	/** Fresnel amplitude transmission coefficient.
	* @param[in] phi angle of incidence in this materials side [rad]
	* @param[in] wn the vacuum wavenumber of the elm wave
	* @param[in] oth the other material from which we refrelct
	* @param[out] T the p polarisation (index 0) and s polarization (index 1) Fresnel coeffs
	*/
	void fresn_coefT(const double phi, const double wn, const DielF & oth, complex<double> T[2]) const
	{
		/* the relative dielectric constant is the ratio of
		* the other materials dielectric constant divided through the one here.
		*/
		const complex<double> nratio = oth.n(wn)/n(wn) ;

		/* Handle the case of perpendicular incidence separately to avoid
		* singularities in the formalism.
		*/
		if (phi == 0.)
			T[0] = T[1] = 2./(nratio+1.) ;
		else
		{
			/* Snell's law in the complex version.
			*/
			const complex<double> psi = asin( sin(phi)/nratio) ;
			T[0] = 2.*sin(psi)*cos(phi)/sin(phi+psi)/cos(phi-psi) ;
			T[1] = 2.*sin(psi)*cos(phi)/sin(phi+psi) ;
		}
	}

	/** Wavenumber range covered by the dielectric functions
	* @param[in,out] w is filled withe the minimum wavenumber in w[0] and the
	*   maximum in w[1] on return, by possibly extending the range provided on input.
	*/
	void wnRange(double w[2])
	{
		for(int i=0; i < tab.size() ; i++)
		{
			w[0] = min(w[0],tab[i].x) ;
			w[1] = max(w[1],tab[i].x) ;
		}
	}

protected:
private:
} ;

/** Report a value to the output stream
* @param[in,out] os the output stream
* @param[in] v the value to be printed
* @return the repositioned output stream
* @author Richard J. Mathar
* @since 2008-05-27
*/
ostream & operator<< (ostream &os, const DielF & v)
{
	for(int i=0 ; i < v.tab.size() ; i++)
		os << v.tab[i] << endl ;
	return os ;
}

/** A layer of material with some thickness and dielectric function.
* @author Richard J. Mathar
* @since 2008-05-27
*/
class Layer : public DielF {
public:
	/* The thickness perpendicular to the surface [m]
	*/
	double thickn ;

	/** Ctor.
	*/
	Layer() : DielF() , thickn(0.)
	{
	}

	/** Ctor.
	* @param[in] diel
	* @param[in] d thickness [m]
	*/
	Layer(const DielF & diel, const double d) : DielF(diel), thickn(d)
	{
	}

	/** Optical path difference by one transit 
	* @param[in] wn vacuum wavenumber
	* @param[in] ang angle of incidence relative to surface normal [rad]
	*   Measured inside the material; no knowledge on refraction at any surface.
	* @return OPD [m], the product of geometric path length and dielectric constant.
	*/
	double opd(const double wn, const double ang) const
	{
		const double nRe = n(wn).real() ;
		return thickn*nRe/cos(ang) ;
	}


protected:
private:
} ;

/** Report a value to the output stream
* @param[in,out] os the output stream
* @param[in] v the value to be printed
* @return the repositioned output stream
* @author Richard J. Mathar
* @since 2008-05-27
*/
ostream & operator<< (ostream &os, const Layer & l)
{
	os << "Layer thickness " << l.thickn << " m\n" ;
	os << "Complex diel. functions and wave numbers:\n" ;
	os << (DielF)l << endl ;
	return os ;
}

/** A beam splitter is a sandwitch of 3 layers: vacumm (air), dielectric surfac, substrate
* @author Richard J. Mathar
* @since 2008-05-27
*/
class Bsplit {
public:
	vector<Layer> combo ;

	Bsplit(const Layer & lint, const Layer & subs)
	{
		/* For the external air interface, we assume good vacuum.
		* Alternatively one would use 1.002 for Paranal standard pressure.
		*/
		const Clist airN(1000, complex<double>(1)) ;
		DielF air ;
		air.disp(airN) ;
		Layer fron(air,1.0) ;
		combo.push_back(fron) ;

		combo.push_back(lint) ;
		combo.push_back(subs) ;
	}

	/** polarization angle between p and s polarization
	* @param[in] R the component R[0] the p polarization complex amplitude, component R[1] the s polarization
	* @return polarization angle [rad]
	*/
	static double pangle(complex<double> R[2])
	{
		return arg(R[1]/R[0]) ;
	}

	/* Amplitude reflection factor for multi-interference reflection, including
	* the phase jumps and optical path differences.
	* @param[in] wn vacuum wavenumber
	* @param[in] ang angle of incidence before the layer
	* @param[in] pol 0 or 1 for p or s polarization
	* @param[out] R amplitude coefficient after reflection from the layer [rad]
	*/
	void Rinterf(const double wn, const double ang, complex<double> R[2], const bool incoh) const
	{
		/* Angle of incidence is ang, angLay is the corresponding angle
		* in the layer, using Snell's law. If h is the thickness of the layer, the sideways motion (geometric) from
		* a passage to the bottom of the layer is h*tan(angLay). Double passage with return gives
		* 2*h*tan(angLay) relative shift at exit. For the right triangle in the directly
		* reflected light this is the hypothenuse, with ang and 90-ang two angles.
		* The OPD for the directly reflected light is the leg x from sin(ang)=x/(2*h*tan()),
		* so x=sin(ang)*2*h*tan(angLay) where sin(angLay)*n'=sin(ang) (Snellius).
		* For the passage (free flight, not yet with the phase jumps at interfaces) we
		* have the geometric path length 2*h/cos(angLay), the OPD 2*h*n'/cos(angLay). Path
		* difference w.r.t. directly refleted (x) is sin(ang)*2*h*tan(angLay)-2*h*n'/cos(angLay)
		* = n'*sin(angLay)*2*h*tan(angLay)-2*h*n'/cos(angLay) = 2*h*n'/cos(angLay)*[sin^2(anglay)-1] = -2*h*n'*cos(anLay)
		* The standard formula of +2*h*n'*cos(angLay) includes already the phase jumps.

		* Evaluation for multi-interference, order N=0,1,2,3...:
		* direct: R12*exp(2*pi*i*wn*x*N) where x=sin(ang)*2*h*tan(angLay)
		* One path through dielectric with d=: 2*h*n'/cos(angLay) the internal OPD, two transmissions,
		*   one internal reflection, and N-1 pieces of the path in the air layer:
		*   T12*exp(2*pi*i*wn*d)*R23*T21*exp{2*pi*wn*x*i*(N-1)}
		* Double path through dielectric
		*   T12*exp(2*pi*i*wn*d)*R23^2*R21*T21*exp(2*pi*wn*d)*exp{2*pi*wn*x*(N-2)}
		* k-fold path through dielectric
		*   T12*exp(2*pi*k*wn*d)*R23^k*R21^(k-1)*T21*exp{2*pi*wn*x*(N-k)}
		* Geometric sum over all complex amplitudes,
		* R12*exp(2*pi*wn*x*N)+T12*T21*exp(2*pi*wn*x*N)/R21*sum_{k=1..infy) exp(2*pi*k*wn*d-2*pi*wn*x*k)*(R23*R21)^k
		* =R12*exp(2*pi*wn*x*N)+T12*T21*exp(2*pi*wn*x*N)/R21*rho/(1-rho)
		* where rho= (R23*R21)*exp(2*pi*i*wn*d-2*pi*i*wn*x); drop the overall phase factor
		* =R12+T12*T21/R21*rho/(1-rho)
		* define delta= 2*pi*wn*(d-x), so rho=R23*R21*exp(i*delta).
		* So far standard Optics see Born or Born&Wolf.
		*/
	
		/* angle of multi-reflection within the dielectric layer
		*/
		const double angLay = combo[0].Tang(ang,wn,combo[1]) ;
	
		/* d and x the optical path differences for the directly reflected and 
		* on the substrate once reflected part [m]. No interface phase jumps applied yet.
		*/
		const double d = 2.* combo[1].opd(wn,angLay) ;
		const double x = 2.* combo[1].thickn* sin(ang)*tan(angLay) ;

		/* phase difference [rad]: factor 100 to convert wavenumber to inverse meters
		* factor 2pi to convert wavenumber to phase
		*/
		const double delta = 200.*M_PI*wn*(d-x) ;

		/* Contribution from the direct reflection on the front surface. We call 1 the air, 2 the
		* coating and 3 the substrate, and the nomenclature has two integers which say where we are
		* before reflecting or transmitting.
		*/
		complex<double> R12[2] ;
		combo[0].fresn_coefR(ang,wn,combo[1],R12) ;

		complex<double> T12[2] ;
		combo[0].fresn_coefT(ang,wn,combo[1],T12) ;

		complex<double> T21[2] ;
		combo[1].fresn_coefT(angLay,wn,combo[0],T21) ;

		complex<double> R23[2] ;
		combo[1].fresn_coefR(angLay,wn,combo[2],R23) ;

		complex<double> R21[2] ;
		combo[1].fresn_coefR(angLay,wn,combo[0],R21) ;

		/* loop over both polarizations
		*/
		for(int pol=0; pol <2 ; pol++)
		{
			const complex<double> rho = R23[pol]*R21[pol]*exp(complex<double>(0.,delta)) ;
	
			/* if incoherent, take only the R12 direct reflection 
			* into account
			*/
			R[pol] = R12[pol] ;
			if ( ! incoh)
				R[pol] += T12[pol]*T21[pol]/R21[pol]*rho/(1.0-rho) ;
		}
	}
	/* Amplitude transmisstion factor for multi-interference reflection, including
	* the phase jumps and optical path differences.
	* @param[in] wn vacuum wavenumber
	* @param[in] ang angle of incidence before the layer
	* @param[out] T phase difference after transmission through layer and substrate [rad]
	*     Component T[0] contains the p, component T[1] the s polarization.
	* @since 2008-05-29
	*/
	void Tinterf(const double wn, const double ang, complex<double> T[2], const bool incoh) const
	{
		/* Angle of incidence is ang, angLay is the corresponding angle
		* in the layer, using Snell's law. If h is the thickness of the layer, the sideways motion (geometric) from
		* a passage to the bottom of the layer is h*tan(angLay). Double passage with return gives
		* 2*h*tan(angLay) relative shift at exit. For the right triangle in the directly
		* reflected light this is the hypothenuse, with ang and 90-ang two angles.
		* The OPD for the directly reflected light is the leg x from sin(ang)=x/(2*h*tan()),
		* so x=sin(ang)*2*h*tan(angLay) where sin(angLay)*n'=sin(ang) (Snellius).
		* For the passage (free flight, not yet with the phase jumps at interfaces) we
		* have the geometric path length 2*h/cos(angLay), the OPD 2*h*n'/cos(angLay). Path
		* difference w.r.t. directly refleted (x) is sin(ang)*2*h*tan(angLay)-2*h*n'/cos(angLay)
		* = n'*sin(angLay)*2*h*tan(angLay)-2*h*n'/cos(angLay) = 2*h*n'/cos(angLay)*[sin^2(anglay)-1] = -2*h*n'*cos(anLay)
		* The standard formula of +2*h*n'*cos(angLay) includes already the phase jumps.

		* Evaluation for multi-interference, order N=0,1,2,3...:
		* Direct through all layers: T12*exp(2*pi*i*wn*d/2)*T23*(const through substr)*T31*exp(2*pi*i*wn*N*x)
		*  where d=2*n(1)*thickness(1)/cos(angLay); first factor contains only d/2 for the single passage through 2.
		* k internal double reflections in the coating:
		* T12*exp(2*pi*i*wn*d/2)*T23*[R23^k R21^k*exp(2*pi*i*wn*d*k)](const through substr)*T31*exp(2*pi*i*wn*(N-k)*x)
		* Geometric sum over k=0..N with D=R23*R21*exp(2*pi*i*wn*d)*exp(-2*pi*i*wn*x):
		* T12*exp(2*pi*i*wn*d/2)*T23*[D^(N+1)-1]/[D-1](const through substr)*T31*exp(2*pi*i*wn*(N-k)*x)
		* In the limit N->infinity, the D^(N+1) vanishes since it contains the product R23*R21 which has a 
		* norm less than 1. Drop constant factors through substrate and exp(2*pi*i*wn*N*x):
		* -T12*T23*T31*exp(2*pi*i*wn*d/2)/[D-1]
		*/
	
		/* angle of multi-reflection within the dielectric layer
		*/
		const double angLay = combo[0].Tang(ang,wn,combo[1]) ;
	
		/* d and x the optical path differences for the directly reflected and 
		* on the substrate once reflected part [m]. No interface phase jumps applied yet.
		*/
		const double d = 2.* combo[1].opd(wn,angLay) ;
		const double x = 2.* combo[1].thickn* sin(ang)*tan(angLay) ;

		/* phase difference [rad]: factor 100 to convert wavenumber to inverse meters
		* factor 2pi to convert wavenumber to phase
		*/
		const double delta = 200.*M_PI*wn*(d-x) ;

		/* Contribution from the direct tranmission on the front surface. We call 1 the air, 2 the
		* coating and 3 the substrate, and the nomenclature has two integers which say where we are
		* before reflecting or transmitting.
		*/
		complex<double> T12[2] ;
		combo[0].fresn_coefT(ang,wn,combo[1],T12) ;

		complex<double> T23[2] ;
		combo[1].fresn_coefT(angLay,wn,combo[2],T23) ;

		/* Angle of passage within the substrate. Instead of transforming angLay once more,
		* we calculate it as if the leight would impinge from the back surface into the substrate.
		*/
		const double angSubs = combo[0].Tang(ang,wn,combo[2]) ;

		complex<double> T31[2] ;
		combo[2].fresn_coefT(angSubs,wn,combo[0],T31) ;

		complex<double> R21[2] ;
		combo[1].fresn_coefR(angLay,wn,combo[0],R21) ;

		complex<double> R23[2] ;
		combo[1].fresn_coefR(angLay,wn,combo[2],R23) ;

		/* loop over both polarizations
		*/
		for(int pol=0; pol <2 ; pol++)
		{
			const complex<double> rho = R23[pol]*R21[pol]*exp(complex<double>(0.,delta)) ;
	
			/* if incoherent, take only the direct transmission into account. delta for the
			* single passage is 2*pi*wn*d/2= pi*wn*d, again with factor 100 to convert Kayser to 1/m.
			*/
			T[pol] = T12[pol]*T23[pol]*T31[pol]*exp(complex<double>(0., 100.*M_PI*wn*d )) ;
			if ( ! incoh)
				T[pol] /= 1.0-rho ;
		}
	}

	/** Wavenumber range covered by the dielectric functions
	* @param[in,out] w is filled withe the minimum wavenumber in w[0] and the
	*   maximum in w[1] on return
	*/
	void wnRange(double w[2])
	{
		/* calls to Layer::wnRange() will extend the range, so we initialize with
		* a much too small upper and a much to large lower bound.
		*/
		w[0]=1.e6 ;
		w[1]= -1. ;
		combo[1].wnRange(w) ;
		combo[2].wnRange(w) ;
	}

protected:
private:
} ;

/** Report a value to the output stream
* @param[in,out] os the output stream
* @param[in] v the value to be printed
* @return the repositioned output stream
* @author Richard J. Mathar
* @since 2008-05-27
*/
ostream & operator<< (ostream &os, const Bsplit & spl)
{
	os << "Beam splitter "  ;
	for(int l=0 ; l < spl.combo.size() ; l++)
	{
		os << "Layer " << l+1 << endl ;
		os << spl.combo[l] << endl ;
	}
	return os ;
}

/** Normalize angle
* @return the angle normalized to the range [-50:310] deg.
* @since 2008-06-09
*/
double normAng(double in)
{
	while ( in < 50.)
		in += 360. ;
	return in ;
}

/** 
* @param[in] prognam
*/
void usage(char *prognam)
{
	cout << "usage: " << prognam << " [-l wavelength/mu | -t diel-layer-thick] [mu1 eps1 eps2 [... mu2 eps1 eps2]] \n";
}

/** Unix interface
* Command line options (must precede the arguments):
*  -l followed by vacuum wavelength of light in micron (default 10 mu)
*  -i followed by an angle [deg] is angle of incidence (default 45 deg)
*  -t followed by reflective surface layer thickness in micron
*        If no value is provided, the program prints a range of values without
*        assuming a specific value.
*  -I with no argument triggers calculation of the incoherent case where
*        only the direct reflection off the first layer is considered,
*        and the multireflecion/interference within the first layer is
*        supressed. (Useful for debugging, Brewster angle etc)
*  -s followed by a positive integer number gives the number of output
*        (sampling) points which are reported. If the -l switch had been
*        used, the sampling is over a range of coating thicknesses at that
*        wavelength, if otherwise the -t had been used, the sampling is
*        over a range of wavelengths keepin the thickness constant.
*        The default is 20.
*
* The dielectric constants (and their dispersion) for the surface layer
* and substrate are given by triplets of numbers. The first triplet
* giving one vacuum wavelength (micron) with two dielectric indices, the
* next triplet giving another vacuum wavelength with another two dielectric
* indices etc. The single wavelength of the -l option takes a Lagrange-interpolation
* between these values to obtain the two dielectric indices of the layer and
* the substrate.
*
* Examples:
*   For vacuum wavelength of 8.9 micron, coating dielectric function 1.1 and
*   substrate dielectric function 1.24 at 9 micron, coating diel. function 1.12
*   and substrate dielectric function 1.26 at 10 micron:
*  MatisseBS -l 8.9 9. 1.1 1.24 10. 1.12 1.26
*
*   For vacuum wavelength of 8.9 micron, surface dielectric function 1.1 and
*   substrate dielectric function 1.24 at 9 micron, surface diel. function 1.12
*   and substrate dielectric function 1.26 at 10 micron, glancing incidence at 20 deg:
*  MatisseBS -l 8.9 -i 20 9. 1.1 1.24 10. 1.12 1.26
*/
int main(int argc, char *argv[])
{

	/* command line option */
	int c;

	/* wavenumber is 10 micron  = 1000 /cm by default */
	double wn = 1.e3 ;

	/* thickness of layer [m] */
	double th = -1. ;

	/* angle of incidence [rad], 45 deg by default */
	double ang = DEG2RAD(45.) ;

	/* a flag to indicate that the multi-interference from the transition
	* from the layer to the substrate is ignored.
	*/
	bool incoh=false ;

	/* a flag to indicate that the ouptut is in pure tabular (uncommented)
	* form, which simplifies plotting
	*/
	bool mute=false ;

	/* Number of sampling points along the wavenumber or thickness range
	*/
	int samp= 20 ;

	while( (c=getopt(argc,argv,"l:t:i:Ims:")) != -1 )
	{
		switch(c)
		{
		case 'l':
			wn= MICRON2WN(atof(optarg)) ;
			break ;
		case 't':
			th= 1.e-6*atof(optarg) ;
			break ;
		case 'i':
			ang = DEG2RAD(atof(optarg)) ;
			break ;
		case 'I':
			incoh = true ;
			break ;
		case 's':
			samp = atoi(optarg) ;
			break ;
		case 'm':
			mute =true ;
			break ;
		case '?':
			usage(argv[0]) ;
			return 1 ;
			break ;
		}
	}

#if 0
	cout << __FILE__ << " " << __LINE__ << endl ;
#endif

	/* dielectric functions of interface and substrate
	*/
	DielF interf ;
	DielF substr ;

	/* collect command line arguments and build dispersion table
	*/
	int dielargs = ( argc-optind)/3 ;
	while (dielargs-- > 0)
	{
		double mu = atof( argv[optind]) ;
		double eps = atof( argv[optind+1]) ;

		Clist tmpI( MICRON2WN(mu), complex<double>(eps,0.) ) ;
		interf.disp(tmpI) ;

		eps = atof( argv[optind+2]) ;
		Clist tmpS( MICRON2WN(mu), complex<double>(eps,0.) ) ;
		substr.disp(tmpS) ;

		optind += 3 ;
	}

	/* parameters of looping over interface thickness
	*/
	double thStrt = (th >0.) ? th : 1.e-7 ;
	double thStep = (th >0.) ? 1. : 1.e-6*WN2MICRON(wn)/40. ;

	Layer back(substr,MATISSE_SUBSTR_THICK) ;

	if ( th < 0 )
	{
		/* parameters of looping over interface thickness
		*/
		double thStop = 5.e-7*WN2MICRON(wn) ;
		double thStrt = (th >0.) ? th : 1.e-7 ;
		double thStep = (th >0.) ? 1. : (thStop-thStrt)/(samp-1) ;

		/* layer thickness not provided by user; build table of reflection
		* coefficients.
		*/
		for(th = thStrt ; th < thStop ; th += thStep )
		{
			Layer top(interf,th) ;
			Bsplit bspl(top,back) ;
	
			/*
			* cout << bspl << endl ;
			*/
	
			/* complex phasors of reflection or transmission; re-cycle R in both cases */
			complex<double> R[2] ;
			bspl.Rinterf(wn,ang,R,incoh) ;
	
			if ( mute) 
			{
				cout << th/1.e-6 << " " << 100.*norm(R[0]) << " " << 100.*norm(R[1]) 
					<< " " << normAng(RAD2DEG(arg(R[0]/R[1]))) ;
			}
			else
			{
			cout << "Coating thickness " << th/1.e-6 << " mu\n" ;
			cout << "Amplitude reflection p-polarization " << R[0] << " " << 100.*norm(R[0]) << " %" << endl ;
			cout << "Amplitude reflection s-polarization " << R[1] << " " << 100.*norm(R[1]) << " %" << endl ;
			cout << "Amplitude reflection p/s ratio " << R[0]/R[1] << " arg " << RAD2DEG(arg(R[0]/R[1])) << " deg" << endl ;
			}
	
			bspl.Tinterf(wn,ang,R,incoh) ;
	
			if ( mute)
			{
				cout << " " << th/1.e-6 << " " << 100.*norm(R[0]) << " " << 100.*norm(R[1]) 
					<< " " << normAng(RAD2DEG(arg(R[0]/R[1]))) << endl ;
			}
			else
			{
			cout << "Amplitude transmission p-polarization " << R[0] << " " << 100.*norm(R[0]) << " %" << endl ;
			cout << "Amplitude transmission s-polarization " << R[1] << " " << 100.*norm(R[1]) << " %" << endl ;
			cout << "Amplitude transmission p/s ratio " << R[0]/R[1] << " arg " << RAD2DEG(arg(R[0]/R[1])) << " deg" << endl ;
			}
		}
	}
	else
	{
		Layer top(interf,th) ;
		Bsplit bspl(top,back) ;

		/* compute the wavenumber range scanned from the dielectric functions known
		*/
		double wrange[2] ;
		bspl.wnRange(wrange) ;
		/* degenerate case: only one wavenumber given; stretch the interval to 2 mu
		*/
		if ( wrange[1] == wrange[0])
		{
			wrange[0] = 1.0/(1.0/wrange[0]+2.e-4) ; 
			wrange[1] = 1.0/(1.0/wrange[0]-2.e-4) ; 
		}

		for( double w = wrange[0] ; w <= wrange[1] ; w += (wrange[1]-wrange[0])/samp )
		{
			/* complex phasors of reflection */
			complex<double> R[2] ;
			bspl.Rinterf(w,ang,R,incoh) ;
	
			if ( mute )
			{
				cout << WN2MICRON(w) << " " << 100.*norm(R[0]) << " " << 100.*norm(R[1]) 
					<< " " << normAng(RAD2DEG(arg(R[0]/R[1]))) ;
			}
			else
			{
			cout << "Wavelength " << WN2MICRON(w) << " mu\n" ;
			cout << "Amplitude reflection p-polarization " << R[0] << " " << 100.*norm(R[0]) << " %" << endl ;
			cout << "Amplitude reflection s-polarization " << R[1] << " " << 100.*norm(R[1]) << " %" << endl ;
			cout << "Amplitude reflection p/s ratio " << R[0]/R[1] << " arg " << RAD2DEG(arg(R[0]/R[1])) << " deg" << endl ;
			}
	
			bspl.Tinterf(wn,ang,R,incoh) ;
	
			if ( mute)
			{
				cout << " " << WN2MICRON(w) << " " << 100.*norm(R[0]) << " " << 100.*norm(R[1]) 
					<< " " << normAng(RAD2DEG(arg(R[0]/R[1]))) << endl ;
			}
			else
			{
			cout << "Amplitude transmission p-polarization " << R[0] << " " << 100.*norm(R[0]) << "%" << endl ;
			cout << "Amplitude transmission s-polarization " << R[1] << " " << 100.*norm(R[1]) << "%" << endl ;
			cout << "Amplitude transmission p/s ratio " << R[0]/R[1] << " arg " << RAD2DEG(arg(R[0]/R[1])) << " deg" << endl ;
			}
		}
	}

	return 0 ;
}