#include <math.h>
#include <stdio.h>
#include <stdlib.h>
/*************************
 Refractive index of dry air, P E Ciddor, Appl Opt 35 (9) (1996)
 usage:
  cidd2 t pair pwat lambdastart lambdaend lambdadelta [xc]
 where t is the temperature in Celsius
 lambdastart, lambdaend and lambdadelta are the wavenlength values (in micrometer)
      to take for evaluation; the output contains values lambdastart<=lambda<lambdaend
 pair: partial pressure of dry air in units of mm Hg
 pwat: partial pressure of water in units of mm Hg

 xc is the CO_2 contents of the air in units of ppm and is set to 350 (ppm) if not given.

 example
  cidd2 12 564.7 2.4 0.2 30. 0.2
 to get Paranal values.

 Richard J. Mathar, http://www.strw.leidenuniv.nl/~mathar 2004-06-16
*****************************/
/*#define GROUP_INDEX */
#define KELVIN 273.15
const double R=8.314510 ;
const double Mw= 0.018015 ;
int main (int argc, char *argv[])
{
	double lambda , lambd[3];	/* wavelength, micrometer */
	double xc ; /* ppm CO2 */
	double k0= 238.0185 ;
	double k1=5792105. ;
	double k2= 57.362 ;
	double k3= 167917 ;
	double Ma ;
	double T,t,p,xw ;
	double Z(double,double,double) ;
	double rho(double,double,double,double,double) ;
	double fenh(double,double) ;
	double Za ;
	double Zw ;
	double Zmoist ;
	double rhoaxs ;
	double rhows ;
	double rhoa , pa ;
	double rhow , pw;

	if( argc < 7)
	{
		printf("usage: %s t p_air/mmHg p_water/mmHg lambda1 lambda2 lambda3 [CO2/ppm]\n",argv[0]) ;
		exit(EXIT_FAILURE) ;
	}
	lambd[0]=atof(argv[4]) ;
	lambd[1]=atof(argv[5]) ;
	lambd[2]=atof(argv[6]) ;
	
	if( argc==8)
		xc=atof(argv[argc-1]) ;
	else
		xc = 350. ;
	Ma= 1.e-3*(28.9635+12.011e-6*(xc-400.)) ;

	pa=atof(argv[2]) ;
	pw=atof(argv[3]) ;

		/* standard air at 15 Celsius, 101325 Pa, zero humidity */
	T=288.15 ;
	p=101325. ;
	xw=0. ;
	Za=Z(p,T,xw) ;

		/* Barell and Sears standard at 20 Celsius, 1333 Pa of water pressure */
	T=293.15 ;
	p=1333 ;
	xw=1. ;
	Zw=Z(p,T,xw) ;

		/* again standard air at 15 Celsius, 101325 Pa, zero humidity */
	T=288.15 ;
	p=101325. ;
	xw=0. ;
	rhoaxs=rho(p,Za,T,xw,Ma) ;

		/* again Barell and Sears standard at 20 Celsius, 1333 Pa of water pressure */
	T=293.15 ;
	p=1333 ;
	xw=1. ;
	rhows=rho(p,Zw,T,xw,Ma) ;

	T=KELVIN+atof(argv[1]) ;
	p= 133.322*(pa+pw); /* paranal: 564.7+2.4 mm Hg */
			/* xw=f*h*svp/p, where
			h*svp is the partial pressure=2.4 mmHg*/
	xw=fenh(p,T)*133.322*pw/p ;
	Zmoist=Z(p,T,xw) ;

	rhoa=p*Ma*(1.-xw)/Zmoist/R/T ;
	rhow=p*Mw*xw/Zmoist/R/T ;

	printf("# Temperature: %e K %e C\n",T,T-KELVIN) ;
	printf("# dry air pressure: %e mmHg %e Pa\n",pa,pa*133.322) ;
	printf("# water pressure: %e mmHg %e Pa\n",pw,pw*133.322) ;
	printf("# total pressure humid air: %e mmHg %e Pa\n",pa+pw,(pa+pw)*133.322) ;
	printf("# Columns from left to right\n") ;
	printf("# 1. wavelength/micron  2. wavenumber/cm-1 3. n-1 dry air STP 4. n-1 humid air 5. n-1 dry air 6. n-1 pure water\n") ;
	for(lambda = lambd[0] ; lambda < lambd[1] ; lambda +=  lambd[2] )
	{
#ifdef GROUP_INDEX
		double ngasx1(double,double) ;
		double ngws1(double) ;
#else
		double naxs1(double,double) ;
		double nws1(double) ;
#endif
#ifdef GROUP_INDEX
		double n8dry=ngasx1(lambda,xc) ;
		double n8wet=ngws1(lambda) ;
#else
		double n8dry=naxs1(lambda,xc) ;
		double n8wet=nws1(lambda) ;
#endif
		double nprop=rhoa/rhoaxs*n8dry+rhow/rhows*n8wet ;
		printf("%e %e %e %e %e %e\n",lambda,1.e4/lambda,n8dry,nprop,rhoa/rhoaxs*n8dry,rhow/rhows*n8wet) ;
	}
	return 0 ;
}

/* apply eq (12) to compute compressibility */
double Z(double p, double T, double xw)
{
	const double t=T-KELVIN ;
	const static double a[3]={1.58123e-6,-2.9331e-8,1.1043e-10} ;
	const static double b[2]={5.707e-6,-2.051e-8} ;
	const static double c[2]={1.9898e-4,-2.376e-6} ;
	const static double d= 1.83e-11 ;
	const static double e= -0.765e-8 ;

	return 1.-p/T*
	(a[0]+a[1]*t+a[2]*t*t+(b[0]+b[1]*t)*xw+(c[0]+c[1]*t)*xw*xw)
		+pow(p/T,2.)*(d+e*xw*xw) ;
}

/* apply eq (4) */
double rho(double p, double thisZ,double T,double xw,double Ma)
{
	return p*Ma/thisZ/R/T*(1.-xw*(1.-Mw/Ma)) ;
}

/* enhancement factor*/
double fenh(double p,double T)
{
	const double t=T-KELVIN ;
	const static double alpha=1.00062 ,
		            beta=3.14e-8 ,
			    gamma=5.6e-7 ;
	return alpha+beta*p+gamma*t*t ;
}

/* equation 1; return nas-1 */
double nas1(double lambda)
{
	double k0= 238.0185 ;
	double k1=5792105. ;
	double k2= 57.362 ;
	double k3= 167917 ;
	double sigma2 = 1./(lambda*lambda) ;	/* sigma^2 */
		/* n8= n-1 in units of 10^8 */
	double n8 = k1/(k0-sigma2) ;

	n8 += k3/(k2-sigma2) ;
	n8 *= 1.e-8 ;
	return n8 ;
}

/* equation 2; return naxs-1 */
double naxs1(double lambda, double xc)
{
	double k0= 238.0185 ;
	double k1=5792105. ;
	double k2= 57.362 ;
	double k3= 167917 ;
	double sigma2 = 1./(lambda*lambda) ;	/* sigma^2 */
		/* n8= n-1 in units of 10^8 */
	double n8 = k1/(k0-sigma2) ;
	n8 += k3/(k2-sigma2) ;
	n8 *= 1.+0.534e-6*(xc-450.) ;
	n8 *= 1.e-8 ;
	return n8 ;
}

/* equation 10; return ngasx-1 */
double ngasx1(double lambda,double xc)
{
	double k0= 238.0185 ;
	double k1=5792105. ;
	double k2= 57.362 ;
	double k3= 167917 ;
	double sigma2 = 1./(lambda*lambda) ;	/* sigma^2 */
		/* n8= n-1 in units of 10^8 */
	double n8 = k1*(k0+sigma2)/pow(k0-sigma2,2.) ;

	n8 += k3*(k2+sigma2)/pow(k2-sigma2,2.) ;
	n8 *= 1.+0.534e-6*(xc-450.) ;
	n8 *= 1.e-8 ;
	return n8 ;
}

/* equation 11, return ngws-1 */
double ngws1(double lambda)
{
	double sigma2 = 1./(lambda*lambda) ;	/* sigma^2 */
			/* note the correction o funits of w[0..3] in
			the Ciddor-Hill followup publication AO 38 p 1663 */
	static const double w[4]={295.235,2.6422,-0.032380,0.004028}  ;
	double n8=1.022*(w[0]+3.*w[1]*sigma2+5.*w[2]*sigma2*sigma2
			+7.*w[3]*pow(sigma2,3.)) ;
	n8 *= 1.e-8 ;
	return n8 ;
}

/* equation 3, return nws-1 */
double nws1(double lambda)
{
	double sigma2 = 1./(lambda*lambda) ;	/* sigma^2 */
	static const double w[4]={295.235,2.6422,-0.032380,0.004028}  ;
	double n8=1.022*(w[0]+w[1]*sigma2+w[2]*sigma2*sigma2 +w[3]*pow(sigma2,3.)) ;
	n8 *= 1.e-8 ;
	return n8 ;
}