// stde.cpp
// Calcularte Permittitity of standard sample
// Parameters are from Tokai Univ.
// Complex functions are AT&T ver. 2.0
// 標準試料の誘電率を計算する
#define STDE_SRC

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <complex.h>
#include "mathconst.h"
#include "stde.h"

/**************************************************************************/
// calculation of Standard Permittivity at given anguler frequency
// w : angular frequency omega
// t : temperature(Celucius degree)
// index:index of array StandardParam StdP[NUM_STD] defined in stde.h
complex<double> STD_Calc(double w, double t, short index)
{
		StandardParam *p;
        complex<double> z;
        static complex<double> es;
        double fmax, lfmax;
        double e00, e0, tau;
        
        p = &StdP[index];		// refer to table in stde.h
        z = complex<double>(0.0, 1.0);
        switch(index){
        	// 真下研と同じ標準計算法
        	case STD_AIR:	/* 分散なし、虚部=0, 実部=1 airの場合を想定 */
        		e0=(double)1;
        		e00=(double)1;
        		tau=(double)0;
        		break;
        	case STD_WATER:
        	case STD_ACETONE:
        	case STD_CHLORO:
        	case STD_BENZENE:
        	case STD_TOLUENE:
        	case STD_ETHANOL:
        	case STD_METHANOL:
        		e00 = p->p[2] * t + p->p[3];
        		//there are two case:frag is 0 or flag is 2
        		if(p->flag == 2){
                	e0 = p->p[0]*t + p->p[1];
                	tau = p->p[4]*t + p->p[5];
        		}
        		else{
                	e0 = 87.740 - 0.40008 * t + 9.398e-4*t*t - 1.410e-6*t*t*t;
                	tau = pow((double)10, (-0.0108*t - 10.790));
        		}
        		break;
        	case STD_WATER_KAATZE:	// Udo Kaatze, J.Chem.Eng.Data 1989, 34, 371-374
        		t += 273.15;		// Kaatze's equation uses temp. in Kelvin.
        		e0 = pow((double)10, 1.94404 - 1.991e-3*(t - 273.15));
				e00 = 5.77 - 2.74e-2 * (t-273.15);
				tau = 3.745e-15 * (1 + 7e-5 * (t-300.65)*(t-300.65)) * exp(2.2957e3 / t);
				break;
			case STD_ETHANOL2:
				t += 273.15;
					// CRC Handbook of Chemistry & Physics, 77th, p6-187
				e0 = 0.15145e3 - 0.87020*t + 0.19570e-2*t*t - 0.15512e-5*t*t*t;  // 163-350(K)
					// NIST, 
				e00 = 4.1871 + 169.27 * exp(-0.028815*t);	// 163-300(K)
					// J.Chem.Phys. vol.104, Feb. 1996, pp.2050
				if(t <= 165){
					lfmax = 8.91 - (double)216/(t-(double)80);	// 108-165(K)
				}
				else{
					lfmax = 12.11 - (double)938/t;	//165-350(K)
				}
				fmax = pow(10.0, lfmax);
				tau = (double)1 / ((double)2 * M_PI * fmax);
				break; 
			case STD_PROPANOL:
				t += 273.15;
					// CRC Handbook of Chemistry & Physics, 77th, p6-188
				e0 = 0.98045e2 - 0.36860*t + 0.36422e-3*t*t;		// 193-350(K)
					// NIST, 
				e00 = 3.0104 + 9.4988*exp(-0.0097533*t);		// 117-350(K)
				if(t <= 136){
					lfmax = 12.4 - (double)806/(t-50.3);		// 107-036(K)
				}
				else if(t >= 182){
					lfmax = 12.89 - (double)1248/t;				// 182-350(K)
				}
				else{
					lfmax = 10.21 - (double)455/(t - (double)73);	// 182-350(K)
				}
				fmax = pow(10.0, lfmax);
				tau = (double)1 / ((double)2 * M_PI * fmax);
				break;
		}
		
		
		// Debye type relaxation
        es = (e0 - e00)/(1.0 + z*w*tau) + e00;
        return es;
}
