/*
# %M% %Y% %I%
# The latest update : %G% at %U%
#
#%Z% llDataCTFCompensate ver %I%
#%Z% Created by 
#%Z%
#%Z% Usage : llDataCTFCompensate
#%Z% Attention
#%Z%
*/
static char __sccs_id[] = "%Z%llDataCTFCompensate ver%I%; Date:%D% %Z%";
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>                  
#define GLOBAL_DECLARATION
#include "../inc/config.h"
#include "llData.h"
#include "String.h"
#include "genUtil.h"
#include "ctfInfo.h"
#define SECTION_MARK " \t,"

extern void llDataCTFSet(llData* llInfo, FILE* fpt);
extern void llDataCTFSet2(llData* llInfo, FILE* fpt);

extern void configFileRead(FILE* fpt, llDataCTFCompensateInfo* info);
extern void configFileUsage(FILE* fpt);

void
main(int argc, char* argv[]) 
{
	long status;
	llDataCTFCompensateInfo info;
    llData   llIn;
    llData   llInfo;

	init0(&info);
    argCheck(&info, argc, argv);
    init1(&info);
	if(info.flagconfigFile) {
		configFileRead(info.fptconfigFile, &info);
	}
    /* Data In */
	llDataFileRead(info.fptIn, &llIn, &status);
    llIn.truePitch = info.truePitch;
    /* llInfo Initialization */
    llInfo = llIn;
    llDataInit(&llInfo);
    lllDataAttributeCopy(&llInfo, &llIn);
    /* CTF Compensation */
	if(info.flagIn1) {
    	llDataCTFSet(&llInfo, info.fptIn1);
	} else if(info.flagIn2) {
		llDataCTFSet2(&llInfo,info.fptIn2);
	} else {
		fprintf(stderr, "You must choose either In1 or In2.\n");
		exit(EXIT_FAILURE);
	}
    lllDataMul(&llIn, &llInfo);
    /* Data Out */
	llDataFileWrite(info.fptOut, &llIn, &status);
	exit(EXIT_SUCCESS);
}

void
additionalUsage()
{
	fprintf(stderr, "----- CTFInformation file 1 -----\n");
	fprintf(stderr, "RR0 F0 PHI0   : If R.lt.RR0, F0=0; If RR0.le.R.lt.RR1, Linear Interpolation F0, F1\n");
	fprintf(stderr, "RR1 F1 PHI1   : If RR1.le.R <RR2,\n"); 
	fprintf(stderr, ".................................................\n");
	fprintf(stderr, "RRn Fn PHIn   : If RRn<= R     , F*0     PHI = 0\n");
	configFileUsage(stderr);
	fprintf(stderr, "----- CTFInformation file 2 -----\n");
	fprintf(stderr, "dF[A] Cs[mm] Acc[kV] ratioOfAmpToPhase\n");
}

void 
llDataCTFSet2(llData* llInfo, FILE* fpt)
{
	ctfInfo ctf;
    llDataParaTypeInteger ll, iR;
    char s[GEN_UTIL_MAX_STRING_LENGTH];
	double CTF, F, R;

	fseek(fpt,0L, SEEK_SET);
    stringGetFromFile(s, "", fpt, NULL, 1);
	ctf.defocus = stringGetNthRealData(s, 1, SECTION_MARK);
	ctf.Cs = stringGetNthRealData(s, 2, SECTION_MARK);
	ctf.kV = stringGetNthRealData(s, 3, SECTION_MARK);
	ctf.ratioOfAmpToPhase = stringGetNthRealData(s, 4, SECTION_MARK);
	ctf.mode = ctfInfoModePhaseAndAmpCTF;

    for(ll=0; ll<llInfo->llMax; ll++) {
        for(iR=0; iR<llInfo->LL[ll].nR; iR++) {
            R = sqrt(SQR(llInfo->LL[ll].data[iR].RR)+SQR(llInfo->LL[ll].nl/llInfo->truePitch));
			CTF = ctfFunction(&ctf, R, 0);
			if(CTF!=0) {
				F = CTF/fabs(CTF);
			} else {
				F = 1.0;
			}
            llInfo->LL[ll].data[iR].F   = fabs(F);
			if(0<F) {
            	llInfo->LL[ll].data[iR].PHI = 0;
			} else {
            	llInfo->LL[ll].data[iR].PHI = M_PI;
			}
            llInfo->LL[ll].data[iR].Re  = F;
            llInfo->LL[ll].data[iR].Im  = 0;
		}
	}
}

void
llDataCTFSet(llData* llInfo, FILE* fpt)
{
    llDataParaTypeInteger ll, iR;
    llDataParaTypeReal    R, RR0, RR1, F0, F1, PHI0, PHI1;
    char s[GEN_UTIL_MAX_STRING_LENGTH];

    fseek(fpt, 0L, SEEK_SET);
    stringGetFromFile(s, "", fpt, NULL, 1);
    RR0  = stringGetNthRealData(s, 1, SECTION_MARK);
    F0   = stringGetNthRealData(s, 2, SECTION_MARK);
    PHI0 = stringGetNthRealData(s, 3, SECTION_MARK)*RADIAN;
    for(ll=0; ll<llInfo->llMax; ll++) {
        for(iR=0; iR<llInfo->LL[ll].nR; iR++) {
            R = sqrt(SQR(llInfo->LL[ll].data[iR].RR)+SQR(llInfo->LL[ll].nl/llInfo->truePitch));
            if(R<RR0) {
                llInfo->LL[ll].data[iR].F   = 0;
                llInfo->LL[ll].data[iR].PHI = 0;
                llInfo->LL[ll].data[iR].Re  = 0;
                llInfo->LL[ll].data[iR].Im  = 0;
            } else {
                break;
            }
        }
    }
    while(NULL!=stringGetFromFile(s, "", fpt, NULL, 1)) {
        RR1  = stringGetNthRealData(s, 1, SECTION_MARK);
        F1   = stringGetNthRealData(s, 2, SECTION_MARK);
        PHI1 = stringGetNthRealData(s, 3, SECTION_MARK)*RADIAN;
        for(ll=0; ll<llInfo->llMax; ll++) {
            for(iR=0; iR<llInfo->LL[ll].nR; iR++) {
                R = sqrt(SQR(llInfo->LL[ll].data[iR].RR)+SQR(llInfo->LL[ll].nl/llInfo->truePitch));
                if(R<RR0) {
                } else if(RR0<=R && R<RR1) {
                    llInfo->LL[ll].data[iR].Re  = (F0*cos(PHI0)*(RR1-R)+F1*cos(PHI1)*(R-RR0))/(RR1-RR0);
                    llInfo->LL[ll].data[iR].Im  = (F0*sin(PHI0)*(RR1-R)+F1*sin(PHI1)*(R-RR0))/(RR1-RR0);
                    llInfo->LL[ll].data[iR].F   = sqrt(SQR(llInfo->LL[ll].data[iR].Re)+SQR(llInfo->LL[ll].data[iR].Im));
                    llInfo->LL[ll].data[iR].PHI = atan2(llInfo->LL[ll].data[iR].Im, llInfo->LL[ll].data[iR].Re);
                } else {
                    break;
                }

            }
        }
        RR0  = RR1;
        F0   = F1;
        PHI0 = PHI1;
    }
    for(ll=0; ll<llInfo->llMax; ll++) {
        for(iR=0; iR<llInfo->LL[ll].nR; iR++) {
            R = sqrt(SQR(llInfo->LL[ll].data[iR].RR)+SQR(llInfo->LL[ll].nl/llInfo->truePitch));
            if(RR0<R) {
                llInfo->LL[ll].data[iR].F   = 0;
                llInfo->LL[ll].data[iR].PHI = 0;
                llInfo->LL[ll].data[iR].Re  = 0;
                llInfo->LL[ll].data[iR].Im  = 0;
            }
        }
    }
}
#include <stdio.h>
#include "String.h"
#include "File.h"
void
configFileRead(FILE* fpt, llDataCTFCompensateInfo* info)
{
    char s[1024];

    fseek(fpt, 0L, SEEK_SET);
    stringGetFromFile(s, "", fpt, stdout, 3);
    if(!info->flagIn) {
        info->In = stringGetNthWord(s, 1, " \t,");
        info->fptIn = fileOpen(info->In, "r");
        info->flagIn++;
    }
    
    stringGetFromFile(s, "", fpt, stdout, 3);
    if(!info->flagIn1) {
        info->In1 = stringGetNthWord(s, 1, " \t,");
        info->fptIn1 = fileOpen(info->In1, "r");
        info->flagIn1++;
    }
    
    stringGetFromFile(s, "", fpt, stdout, 3);
    if(!info->flagOut) {
        info->Out = stringGetNthWord(s, 1, " \t,");
        info->fptOut = fileOpen(info->Out, "w");
        info->flagOut++;
    }
    
    stringGetFromFile(s, "", fpt, stdout, 3);
    if(!info->flagtruePitch) {
        info->truePitch = stringGetNthRealData(s, 1, " \t,");
        info->flagtruePitch++;
    }
    
    stringGetFromFile(s, "", fpt, stdout, 3);
    if(!info->flagmode) {
        info->mode = stringGetNthIntegerData(s, 1, " \t,");
        info->flagmode++;
    }
    
}
#include <stdio.h>
#include "String.h"
#include "File.h"
void
configFileUsage(FILE* fpt)
{
    fprintf(fpt, "----- Control File -----\n");
    fprintf(fpt, "In[inFile:NULL]\t\t\t\t#InputDataFile\n");
    fprintf(fpt, "In1[inFile:NULL]\t\t\t\t#InputCTFInformation\n");
    fprintf(fpt, "Out[outFile:NULL]\t\t\t\t#OutputDataFile\n");
    fprintf(fpt, "truePitch[Real:100.0]\t\t\t\t#truePitch[A]\n");
    fprintf(fpt, "mode[Integer:0]\t\t\t\t#Mode\n");
}
