/** read filter CAL data **/
/** for reading the format of 2010 data **/
/** 2011.4/27 by Akira Kadokura **/

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

static char in_dir[] ="/data/aurora/amuaurora/Syowa/Syowa_WMI/CAL/Syowa_WMI_filter_CAL_20221003";


main(argc,argv)
    int argc;
    char *argv[];
{
	FILE *fin,*fout,*fopen();
	char fin_name[100],fout_name[100];
	static char ck[9],cmmdd[5],cmm[3];

	char chead[256],x[10],T[10];
	int i,iend,icenter;
	float hacho[2000],touka[2000],EW,EW2,EW20;
	float rhacho,rtouka,chacho,ctouka,ref,hacho1,hacho2,FWHM;

  /******************/
    if(argc != 3){
		printf("Please key-in (e.g. rdfilterCAL (file name) (wave length to calculate T)\n");
        exit(0);
    };

	sprintf(fin_name,"%s/%s",in_dir,argv[1]);
    fin = fopen(fin_name,"r");
	if(fin==0){
		printf("%s does not exist !!\n",fin_name);
		exit(0);
	};

	sprintf(fout_name,"%s/%s.smry.txt",in_dir,argv[1]);
    fout = fopen(fout_name,"w");

  /** dummy read **/
  	while(1){
		fgets(chead,255,fin);
		if(chead[0]=='n' && chead[1]=='m') break;
	};

  /** read wave_length & transmission rate data **/
  	EW = 0.0; /** Equivalent Width **/

	i=-1;
	while(fgets(chead,255,fin)!=NULL){
		++i;
		sscanf(chead,"%s %s",x,T);
		hacho[i] = atof(x);
		touka[i] = atof(T)/100.0;
		if(i!=0) EW += (hacho[i-1]-hacho[i])*(touka[i]+touka[i-1])/2.0;
/*
		printf("%7.3f %8.6f\n",hacho[i],touka[i]);
*/
    };
	iend = i;


  /** calculate T at the assigned wave length **/
	rhacho = atof(argv[2]);
	for(i=1;i<=iend-1;i++){
		if(rhacho<=hacho[i] && rhacho>hacho[i+1]){
			rtouka = (rhacho - hacho[i+1])/(hacho[i] - hacho[i+1])*(touka[i] - touka[i+1]) + touka[i+1];
		    break;
		};
	};

  /** calculate center wave length **/
  	EW2 = 0.0;
	for(i=1;i<=iend;i++){
		EW2 += (hacho[i-1]-hacho[i])*(touka[i]+touka[i-1])/2.0;
		if(EW2>=EW/2.0){
			icenter = i;
			chacho = (EW20 - EW/2.0)/(EW20 - EW2)*(hacho[i-1] - hacho[i]) + hacho[i];
			ctouka = (EW20 - EW/2.0)/(EW20 - EW2)*(touka[i-1] - touka[i]) + touka[i];
			break;
		};
		EW20 = EW2;
	};

  /** calculate FWHM **/
  	ref = ctouka/2.0;
	for(i=0;i<=iend-1;i++){
		if(touka[i]<ref && touka[i+1]>=ref){
			hacho1 = (ref - touka[i+1])/(touka[i]-touka[i+1])*(hacho[i]-hacho[i+1]) + hacho[i+1];
		};
		if(touka[i]>=ref && touka[i+1]<ref){
			hacho2 = (ref - touka[i+1])/(touka[i]-touka[i+1])*(hacho[i]-hacho[i+1]) + hacho[i+1];
			FWHM = fabs((double)(hacho2 - hacho1));
			if(FWHM<0.1){
				continue;
			}else{
				break;
			};
		};
	};

	printf("hacho1:%8.3f hacho2:%8.3f FWHM:%7.4f\n",hacho1,hacho2,FWHM);
	printf("center(nm):%8.3f T(center):%7.4f FWHM(nm):%7.4f EW(nm):%7.4f\n",chacho,ctouka,FWHM,EW);
	printf("rhacho(nm):%8.3f rtouka:%7.3f\n",rhacho,rtouka);
	printf("\n");

	fprintf(fout,"center(nm):%8.3f T(center):%7.4f FWHM(nm):%7.4f EW(nm):%7.4f\n",chacho,ctouka,FWHM,EW);

	fclose(fin);
	fclose(fout);
};