// TDR_calc.cp
// Calculation routines for TDR spectroscopy
// TDRの反射波のデータのフーリエ変換とイテレーションを行う
//            by apj
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <Strings.h>
#include <complex.h>

#include "debugs.h"
#include "inputdialog.h"
#include "datautils.h"
#include "TDR.h"
#include "stde.h"

// global for calcdialog
char curr_fname[64];
double curr_freq;
Boolean cancelled;

// 差スペクトルを微分するかどうか。従来法でやるときはこれをコメントにする。
//#define USE_DIFFERENTIATE
//#define CHECK    //rhos出力スイッチ

// Top routine of calculation
// Plase call this routine after all data is read and all parameter is given
// Use global variable in TDR.h
short StartCalc(TDR_PARAMETER *pr)
{
	short i;
	double baseav, tailav, height, dbh, autodc;   // for conductivity calculation
	double w;
	complex<double> z, ft_dR, ft_Rdt; 
	complex<double> ze;
	complex<double> epss, epsx;
	complex<double> rho;
	double sigbe; 	/* sigma/eps0 */
	short nit;
	DialogPtr dlg;
	short itemtype, ans;
	Handle itemH;
	char strlf[256];
	GrafPtr cptr;
	Rect box;
	Str255 fnPstr;
#ifdef CHECK
	FILE *fp;
#endif

	// Do Calibration?
	if((pr->baseCs == 0 && pr->baseCe == 0) || (pr->cutCs == 0 && pr->cutCe == 0) ){
		pr->calibflag = false;
	}
	else{
		pr->calibflag = true;
	}

	if(pr->calibflag == true){	// DC Calib.　微分で計算するときには無意味。NaClを使う場合
		Data_difference(RB, RC, dR, N_DATA);
		Data_differentiate(TA, RB, Rdt, N_DATA);

		if(SetBaseLine(dR, pr->baseCs, pr->baseCe, N_DATA) != 0) return -1;
		if(SetBaseLine(Rdt,  pr->baseCs, pr->baseCe, N_DATA) != 0) return -1;
	
		if(CutData(dR, pr->cutCs, pr->cutCe, N_DATA) != 0) return -1;
		if(CutDataSimple(Rdt, pr->cutCs, pr->cutCe, N_DATA) != 0) return -1;
		
		pr->n_cdata = pr->cutCe - pr->baseCs;
	
		pr->dcc = *(dR + pr->cutCe);		// last data of cutting keeps difference from zero_line
		Data_add(RB, RC, Rpp, N_DATA);
		if(Data_average(Rpp, &baseav, pr->baseCs, pr->baseCe, N_DATA) != 0) return -1;
		if(Data_average(Rpp, &tailav, pr->cutCs, pr->cutCe, N_DATA) != 0) return -1;
		height = tailav - baseav;
		dbh = pr->dcc / height;
	
		// Causion! Following routines are varid only the case of no calibration.
		// when you do manual calibration, you must input G0 value to calib.conductibity
		pr->gC = 0.0;
		autodc = (pr->gC + dbh) / (1.0 + pr->gC * dbh) * Cv / pr->gd;
		pr->gC  = dbh;				/* conductivity */
		pr->sbec = autodc;
	}
	

	Data_difference(RS, RX, dR, N_DATA);
#ifdef USE_DIFFERENTIATE
	Data_differentiate(TA, dR, dR_dt, N_DATA);
#endif
	Data_differentiate(TA, RS, Rdt, N_DATA);

	if(SetBaseLine(dR, pr->baseXs, pr->baseXe, N_DATA) != 0) return -1;
	if(SetBaseLine(Rdt,  pr->baseXs, pr->baseXe, N_DATA) != 0) return -1;
	
	if(CutData(dR, pr->cutXs, pr->cutXe, N_DATA) != 0) return -1;
	if(CutDataSimple(Rdt, pr->cutXs, pr->cutXe, N_DATA) != 0) return -1;

#ifdef USE_DIFFERENTIATE
	if(SetBaseLine(dR_dt, pr->baseXs, pr->baseXe, N_DATA) != 0) return -1;
	if(CutData(dR_dt, pr->cutXs, pr->cutXe, N_DATA) != 0) return -1;
#endif

#ifdef CHECK
	fp = fopen("diffdata", "wt");
	for(i = 0 ; i < N_DATA ; i++){
		fprintf(fp, "%e\n", *(dR_dt + i));
	}
	fclose(fp);
#endif CHECK


	pr->n_data = pr->cutXe - pr->baseXs;
	
	pr->dcx = *(dR + pr->cutXe);		// last data of cutting keeps difference from zero_line
	Data_add(RS, RX, Rpp, N_DATA);
	if(Data_average(Rpp, &baseav, pr->baseXs, pr->baseXe, N_DATA) != 0) return -1;
	if(Data_average(Rpp, &tailav, pr->cutXs, pr->cutXe, N_DATA) != 0) return -1;
	height = tailav - baseav;
	dbh = pr->dcx / height;
	
	// Causion! Following routines are good only the case of no calibration.
	// when you do manual calibration, you must input G0 value to calib.conductibity
	autodc = (pr->gC + dbh) / (1.0 + pr->gC * dbh) * Cv / pr->gd;
	pr->gX  = dbh;				/* conductivity */
	pr->sbex = autodc;
	
	
	// make frequency axis for F.T.
	// check frequency range
	if(pr->lfe < pr->lfs){ 
		TDR_paramerr("\pEnd frequency is larger than Start frequency.", "\p","\p","\p");
		return -1;
	}
	if(pr->dlf <= 0){
		TDR_paramerr("\pStep frequency must be larger than zero.", "\p","\p","\p");
		return -1;
	}
	
	i = 0;
	while(1){
		*(lfa + i) = pr->lfs + (double)i * pr->dlf;
		*(fa + i) = pow(10.0, *(lfa + i));
		if(*(lfa + i) > pr->lfe) break;
		i++;
		if(i >= N_DATA){
			ans = CautionAlert(TDR_ALERT_TOOMANY_FREQ, 0);
			if(pr->launch_flag == LAUNCH_BY_DD) return -1;
			if(ans == 1){ /* cancelled */
				return -1;
			}
		}

	}
	pr->n_spec = i;

	GetPort(&cptr);
	cancelled = false;
	dlg = SetupCalcDlg();
	// set current data file name to dialog
	GetDialogItem(dlg, TDR_CALC_DITEM_FILE, &itemtype, &itemH, &box);
	strncpy((char *)fnPstr, (char *)(current_fss->name), 256);
	SetDialogItemText(itemH, fnPstr);
	DrawDialog(dlg);

#ifdef CHECK
	fp=fopen("rhotest", "wt");
#endif

	for(i = 0 ; i < pr->n_spec ; i++){
		// show message
		GetDialogItem(dlg, TDR_CALC_DITEM_FREQ, &itemtype, &itemH, &box);
		sprintf(strlf, "%5.3le", *(lfa + i));
		SetDialogItemText(itemH, C2PStr(strlf));
		// cancell?
		if(cancelled == true){
			SetPort(cptr);
			return -1;
		}		
		
		// Fourier Transform
		w = (double)2 * M_PI * (*(fa + i));
//		ft_dR = ComplexFT(TA, dR, pr->baseXs, pr->cutXe, pr->dt, w);
//		ft_Rdt= ComplexFT(TA, Rdt, pr->baseXs, pr->cutXe, pr->dt, w);
		ft_dR = ComplexFT(TA, dR, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
		ft_Rdt= ComplexFT(TA, Rdt, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
#ifdef USE_DIFFERENTIATE
		ft_dR = ComplexFT(TA, dR_dt, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
		z = complex<double>(0.0, 1.0);
		ft_dR = (double)-1 * ft_dR * z / w / (*(TA+1) - *(TA));
#endif		

#ifdef CHECK
		fprintf(fp, "%e %e\n", real(ft_dR), imag(ft_dR));
#endif
		// Add DC
		z = complex<double>(0.0, 1.0);
		ze = exp(-z * w * (*(TA + pr->cutXe)));
		ft_dR += (pr->dcx / (z * w)) * ze;
		ft_Rdt *= exp(z * w * (*(TA+1) - *(TA)) / (double)2.0) / ((double)2 * z * sin(w * (*(TA+1) - *(TA))/(double)2));
		ft_Rdt = (double)2 * ft_Rdt - ft_dR;
		
		// get standard permittivity
		epss = STD_Calc(w, pr->temperature, pr->id_standard);
		*(standard_eps + i) = epss;			/* for debug */
		epss += (pr->gC * Cv) / (z * w * pr->gd);
		
		// Iteration
		rho = ft_dR / ft_Rdt;
		epsx = Do_Iteration(rho, w, pr, epss, &nit);
				
		// calcurate sample's permittivity
		epsx = epsx + z * pr->sbex / w;
		*(spect + i) = complex<double>(real(epsx), -imag(epsx));
		*(n_it + i) = nit;
		
	}
#ifdef CHECK
	fclose(fp);
#endif

	
	DisposeCalcDlg(dlg);
	
	SetPort(cptr);

	// automatic calibraton for low frequency
	if(pr->DCf < 0.1){		/* no calibration */
		for(i = 0 ; i < pr->n_spec ; i++) *(calimag + i) = 0.0;		/* clear */
		return 0;
	}

	// Fourier Transform
	w = (double)2 * M_PI * pr->DCf;
	ft_dR = ComplexFT(TA, dR, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
#ifdef USE_DIFFERENTIATE
	ft_dR = ComplexFT(TA, dR_dt, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
	z = complex<double>(0.0, 1.0);
	ft_dR = (double)-1 * ft_dR * z / w / (*(TA+1) - *(TA));
#endif		
	ft_Rdt= ComplexFT(TA, Rdt, pr->baseXs, pr->cutXe, (*(TA+1) - *(TA)), w);
		
	// Add DC
	z = complex<double>(0.0, 1.0);
	ze = exp(-z * w * TA[pr->cutXe]);
	ft_dR += (pr->dcx / (z * w)) * ze;
	ft_Rdt *= exp(z * w * (*(TA+1) - *(TA)) / (double)2.0) / ((double)2 * z * sin(w * (*(TA+1) - *(TA))/(double)2));
	ft_Rdt = (double)2 * ft_Rdt - ft_dR;
		
	// get standard permittivity
	epss = STD_Calc(w, pr->temperature, pr->id_standard);
	*(standard_eps + i) = epss;			/* for debug */
	epss += (pr->gC * Cv) / (z * w * pr->gd);
		
	// Iteration
	rho = ft_dR / ft_Rdt;
	epsx = Do_Iteration(rho, w, pr, epss, &nit);
				
	// calcurate sample's permittivity
	epsx = epsx + z * pr->sbex / w;

	sigbe = (double)-1 * imag(epsx) * w;
	for(i = 0 ; i < pr->n_spec ; i++){
		w = (double)2 * M_PI * (*(fa + i));
		*(calimag + i) = imag(*(spect + i)) - sigbe/w;
	}

	return 0;
}

	

// calculate x1-x2 and return result to x
// the number of data points is n
short Data_difference(double *x1, double *x2, double *x, short n)
{
	short i;
	
	for(i = 0 ; i < n ; i++){
		*(x + i)  = *(x1 + i) - *(x2 + i);
	}
	return 0;
}


// differential calculation
// dx is (dy/dx)
// the number of data points is n
short Data_differentiate(double *x, double *y, double *dy, short n)
{
	short i;
	
	for(i = 2 ; i < n ; i++){
		*(dy + i) = *(y + i) - *(y + i - 1);
	}
	return 0;
}

// calculate x1 + x2 and return result to x
// the number of data points is n
short Data_add(double *x1, double *x2, double *x, short n)
{
	short i;
	
	for(i = 0 ; i < n ; i++){
		*(x + i) = *(x1 + i) + *(x2 + i);
	}
	return 0;
}

// calculate average value of data
// averaging from start to end
// return average in *av
// the number of data is n
short Data_average(double *x, double *av, short start, short end, short n)
{
	short i, dn;
	
	if(start >= end || start < 0 || end >= n){
		TDR_paramerr("\pIllegal start point and end point","\pCannot calcurate average","\p","\p");
		return -1;
	}
	*av = 0.0;
	dn = end - start + 1;
	
	for(i = start ; i <= end ; i++){
		*av += *(x + i);
	}
	*av = *av / (double)dn;
	return 0;
}

// Set baseline to zero
// average of data from start to end is baseline.
// *data changed because of subtract baseline
// the number of data points is n
short SetBaseLine(double *data, short start, short end, short n)
{
	short i, dnum;
	double av = 0.0;
	
	if (end <= start || start < 0 || end > n) {
		TDR_paramerr("\pIllegal start point and end point","\pCannot calcurate baseline.", "\p","\p");
		return -1;
	}
	
	// calculate average
	dnum = end - start + 1;
	for(i = start ; i <= end ; i++){
		av += *(data + i);
	}
	av /= (double)dnum;
	
	// set filter of startpoint
	for(i = 0 ; i < start ; i++){
		*(data + i) = 0.0;
	}
	for(i = start ; i < end ; i++){
		*(data + i) = (*(data + i) - av) * (double)(i - start) / (double)dnum;
	}
	for(i = end ; i < n ; i++){
		*(data + i) -= av;
	}
	return 0;
}


// cutting data
// the number of data points is n
// the last data of array *data keeps 'DC'
short CutData(double *data, short start, short end, short n)
{
	short i, dnum;
	double av = 0.0;
	
	
	if (end <= start || start < 0 || end > n){
		TDR_paramerr("\pIllegal start point and end point","\pCannot cut data.","\p","\p");
		return -1;
	}
	// calculate average
	dnum = end - start + 1;
	for(i = start ; i <= end ; i++){
		av += *(data + i);
	}
	av /= (double)dnum;

	for(i = start ; i < end ; i++){
		*(data + i) = (*(data + i) - av) * (double)(end - i) / (double)dnum + av;
	}
	for(i = end ; i < n ; i++){
		*(data + i) = av;
	}
	return 0;
}

// cutting data (simple case)
// the number of data points is n
// the last data of array *data is 0.0
short CutDataSimple(double *data, short start, short end, short n)
{
	short i, dnum;
	double av = 0.0;
	
	if (end <= start || start < 0 || end > n){
		TDR_paramerr("\pIllegal start point and end point","\pCannot cut data.","\p","\p");
		return -1;
	}
	
	dnum = end - start + 1;

	for(i = start ; i < end ; i++){
		*(data + i) = *(data + i) * (double)(end - i) / (double)dnum;
	}
	for(i = end ; i < n ; i++){
		*(data + i) = av;
	}
	return 0;
}

// Complex Fourier Transform
// *t ; time axis, *x : data
// start, end, range of transform
// dt : time interval
// w  : anguler frequency of F.T.
// return commplex FT result.
complex<double> ComplexFT(double *t, double *x, short start, short end, double dt, double w)
{
	static complex<double> result;
	short i;
	complex<double> z = complex<double>(0.0, 1.0);
	double ftreal, ftimag;
	
	ftreal = 0.0;
	ftimag = 0.0;
	for(i = start ; i <= end; i++){
		ftreal += *(x + i) * cos(w * (*(t + i)));
		ftimag += *(x + i) * sin(w * (*(t + i)));
	}
	result = complex<double>(ftreal, -ftimag);	// exp(-iwt)を掛けている。また、２πで規格化していない。
	result *= dt;
	
	return result;
}

// Iteration
// return sample's permittivity
// and return the number of times of ireration in n_it
complex<double> Do_Iteration(complex<double> rho, double w, TDR_PARAMETER *pr, complex<double> es, short *n_it)
{
	static complex<double> ex;
	complex<double> z = complex<double>(0.0, 1.0);
	complex<double> zw, fs, fa, esg, er, f, ff, ef;
	double rzf;
	
	zw = z * w;
	fs = ZcotZ(w, pr->d, es);
	fa = Cv * fs / (zw * pr->gd * es);
	
	ef = es * (1.0 + rho*fa)/(1.0+rho/fa);
	f = ZcotZ(w, pr->d, ef)/fs;
	ex = ef * f;
	rzf = abs(f-1.0);
	
	*n_it = 0;
	
	while(rzf >= pr->condition_done){
		ff = f;
		f = ZcotZ(w, pr->d, ex) / fs;
		ex = ef * f;
		rzf = abs(f - ff);
		*n_it += 1;
		if(*n_it > pr->limit_nit) break;
	}
	return ex;
}


// function ZcotZ
complex<double> ZcotZ(double w, double d, complex<double> es)
{
	static complex<double> result;
	complex<double> z = complex<double>(0.0, 1.0);
	complex<double> x, x2;
	double rx;
	
	x = w * d / Cv * sqrt(es);
	x2 = w * w * d * d / Cv / Cv * es;
	rx = abs(x);
	
	if(rx < 1.0e-1){
		result = 1.0 - x2 /  3.0 - x2 * x2 / 45.0;
	}
	else{
		result = (x / sin(x)) * cos(x);
	}
	return result;
}

/********************************************/
short CalcStandard(TDR_PARAMETER *pr)
{
	int i;
	int ans;
	double w;
	
	i = 0;
	while(1){
		*(lfa + i) = pr->lfs + (double)i * pr->dlf;
		*(fa + i) = pow(10.0, *(lfa + i));
		if(*(lfa + i) > pr->lfe) break;
		i++;
		if(i >= N_DATA){
			ans = CautionAlert(TDR_ALERT_TOOMANY_FREQ, 0);
			if(ans == 1){ /* cancelled */
				return -1;
			}
		}
	}
	pr->n_spec = i;

	for(i = 0 ; i < pr->n_spec ; i++){
		w = (double)2 * M_PI * (*(fa + i));
		*(standard_eps + i) = STD_Calc(w, pr->temperature, pr->id_standard);

	}
	return 0;
}