#pragma hdrstop
#pragma argsused

#ifdef _WIN32
#include <tchar.h>
#else
  typedef char _TCHAR;
  #define _tmain main
#endif

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <conio.h>

double AbsCoef(double nu, double T0, double F)
{
	nu=nu*1e12;
	const double cc0=0.002423;
	const double cc1T0=-5.714e-08;
	const double cc2T0=5.037e-10;
	const double cc1nu=2.263e+06;
	const double c0=2.99792458e8;
	const double a[]={79.23882, 3.815866, 1.634967};
	const double b[]={0.004300598, 0.01117295, 0.006841548};
	const double c[]={1.382264E-13, 3.510354E-16, 6.30035E-15};
	const double d[]={652.7648, 1249.533, 405.5169};
	const double Tc=133.1383;
	const double p[]={0.8379692, -0.006118594, -0.000012936798, 4235901000000.0,
	-14260880000.0, 273815700, -1246943, 9.618642E-14, 1.795786E-16, -9.310017E-18,
	1.655473E-19, 0.6165532, 0.007238532, -0.00009523366, 15983170000000.0,
	-74413570000.0,497448000, 2.882476E-14, -3.142118E-16, 3.528051E-18};
	double es=87.9144-0.404399*T0+9.58726E-4*T0*T0-1.32802e-6*T0*T0*T0;
	double Delta_1=a[0]*exp(-b[0]*T0);
	double Delta_2=a[1]*exp(-b[1]*T0);
	double Delta_3=a[2]*exp(-b[2]*T0);
	double tau_1=c[0]*exp(d[0]/(T0+Tc));
	double tau_2=c[1]*exp(d[1]/(T0+Tc));
	double tau_3=c[2]*exp(d[2]/(T0+Tc));
	double Delta_4=p[0]+p[1]*T0+p[2]*T0*T0;
	double f_1=p[3]+p[4]*T0+p[5]*T0*T0+p[6]*T0*T0*T0;
	double tau_4=p[7]+p[8]*T0+p[9]*T0*T0+p[10]*T0*T0*T0;
	double Delta_5=p[11]+p[12]*T0+p[13]*T0*T0;
	double f_2=p[14]+p[15]*T0+p[16]*T0*T0;
	double tau_5=p[17]+p[18]*T0+p[19]*T0*T0;
	double K1_1=f_1+nu;
	double K1_2=f_1-nu;
	double K1_3=f_2+nu;
	double K1_4=f_2-nu;
	double C1_1=pow(tau_1,2)*Delta_1/(1+pow(2*M_PI*nu*tau_1,2));
	double C1_2=pow(tau_2,2)*Delta_2/(1+pow(2*M_PI*nu*tau_2,2));
	double C1_3=pow(tau_3,2)*Delta_3/(1+pow(2*M_PI*nu*tau_3,2));
	double C2_1=nu*K1_1/(1+pow(2*M_PI*tau_4*K1_1,2));
	double C2_2=nu*K1_2/(1+pow(2*M_PI*tau_4*K1_2,2));
	double C2_3=nu*K1_3/(1+pow(2*M_PI*tau_5*K1_3,2));
	double C2_4=nu*K1_4/(1+pow(2*M_PI*tau_5*K1_4,2));
	double C3_1=pow(2*M_PI*tau_4,2)*Delta_4/2;
	double C3_2=pow(2*M_PI*tau_5,2)*Delta_5/2;
	double C4_1=tau_1*Delta_1/(1+pow(2*M_PI*nu*tau_1,2));
	double C4_2=tau_2*Delta_2/(1+pow(2*M_PI*nu*tau_2,2));
	double C4_3=tau_3*Delta_3/(1+pow(2*M_PI*nu*tau_3,2));
	double C5_1=1/(1+pow(2*M_PI*tau_4*K1_1,2));
	double C5_2=1/(1+pow(2*M_PI*tau_4*K1_2,2));
	double C5_3=1/(1+pow(2*M_PI*tau_5*K1_3,2));
	double C5_4=1/(1+pow(2*M_PI*tau_5*K1_4,2));
	double C6_1=M_PI*nu*tau_4*Delta_4;
	double C6_2=M_PI*nu*tau_5*Delta_5;
	double C0=2*M_PI*nu;
	double e1=es-pow(C0,2)*(C1_1+C1_2+C1_3)-C3_1*(C2_1-C2_2)-C3_2*(C2_3-C2_4);
	double e2=C0*(C4_1+C4_2+C4_3)+C6_1*(C5_1+C5_2)+C6_2*(C5_3+C5_4);
	return (4*M_PI*nu/c0*sqrt(sqrt(e1*e1+e2*e2)/2-e1/2))*
					 exp(((-4)*M_PI*nu/c0*sqrt(sqrt(e1*e1+e2*e2)/2-e1/2))/cc1nu)*
					 (cc0+cc1T0*T0+cc2T0*T0*T0)*F;
}

#define PAUSE 0
#define TABLE 1
#define HNDTB 2

 int _tmain(int argc, _TCHAR* argv[])
{
	#define NUMITEMS 3
	char *inputss[] = {"Frequency, THz: ", "Temperature, øC: ", "Fluence, J/cm^2: "};
	double vals[NUMITEMS], res;
	char string[80], *str, *pstr;
	bool keys[]={true, false, false};
newcalc:
	if (argc<NUMITEMS){
		for(int i=0;i<NUMITEMS; i++)
		{
			while (1)
			{
				printf("%s",inputss[i]);
				str=gets_s(string,20);
				if(str==NULL)
				{
					printf("Error read value");
					return -2;
				}
				if(strcmpi(str,"exit")==0)
					return -1;
				try{vals[i]=strtod(str,&pstr);}
				catch(...)
				{
					printf("\nError convert value <%s> to number.\n",str);
					printf("\nPress any key for exit...\n");
					while (!kbhit()){}
					int c = getch();
					return -2;
				}
				if(*str!=*pstr)break;
				printf("\nValue must be a number!\nOr enter 'exit'\n\n");
			}
		}
	}else{
		for(int i=0; i < NUMITEMS; i++)
		{
			try
			{
				vals[i]=strtod(argv[i+1],&pstr);
				if (*argv[i+1]==*pstr)
				{
					printf("\nError parameter %d. Value <%s> must be a number.\n",i,argv[i+1]);
					printf("\nPress any key for exit...\n");
					while (!kbhit()){}
					int c = getch();
					return -2;
				}
			}
			catch(...)
			{
				printf("Error parameter %d, value <%s>.\n",i,argv[i+1]);
				printf("\nPress any key for exit...\n");
				while (!kbhit()){}
				int c = getch();
				return -2;
			}
		}
		if(argc>(NUMITEMS+1)) {
			for(int i=NUMITEMS+1; i<argc; i++)
			{
				switch (argv[i][0])
				{
					case '-':
						if(strcmpi(argv[i],"--p")==0)
						{
							keys[PAUSE]=false;
						}else if(strcmpi(argv[i],"-t")==0){
							keys[TABLE]=true;
						}else if(strcmpi(argv[i],"-h")==0){
							keys[TABLE]=true;
							keys[HNDTB]=true;
						}
					break;
					default:;
				}
			}
		}
	}

	res=AbsCoef(vals[0], vals[1], vals[2]);
	if(keys[HNDTB])
	{
		printf("Freq.,THz|Temp.0,øC| Fl.,J/cm^2|    dT\n", vals[0], vals[1], vals[2], res);
		printf("---------+---------+-----------+----------\n");
		printf(" % -7.2f | % -6.2f  | % -7.4g | % -8.4g\n", vals[0], vals[1], vals[2], res);
	}
	else if(keys[TABLE])
	{
		printf(" % -7.2f | % -6.2f  | % -7.4g | % -8.4g\n", vals[0], vals[1], vals[2], res);
	}
	else
	{
		printf("\nt = %g\n",res);
	}

	if(keys[PAUSE])
	{
		printf("\nPress any key for exit or 'Enter' for recalculate...\n\n");
		while (!kbhit()){}
		int c = getch();
		if (c==13)
		{
			argc=0;
			keys[HNDTB]=false;
			goto newcalc;
		}
	}
	return 0;
}