/*
Information
Choose the directory indicated when prompted.

Again this is a text image which may be easily opened in other software.
*/

var calculateUncorrected=false;

macro"PS_homoFRET_r"{
run("Set Measurements...", "redirect=None decimal=9");

dir0=getDirectory("Pick the directory parallel and perpendicular data")
dir0List=getFileList(dir0);
for(i=0;i<dir0List.length;i++){
	if(findStringFragment(dir0List[i],"_Parallel_raw_values_"))
	ParallelOpenFilePath=dir0+dir0List[i];
	if(findStringFragment(dir0List[i],"_Perpendicular_raw_values_"))
	PerpendicularOpenFilePath=dir0+dir0List[i];
}

gFactor=1.054;
NAObj=1.4;//objective
indexRefraction=1.52;//assume oil
R0=5.2;//Dronpa-Dronpa R0
fitFractionCutOff=0.8;//to avoid the noisy data at the end of the anisotropy curves; it will use from 0 to (1-fitFraction) for fitting


//******************** put the Parallel data into arrays 
run("Text Image... ", "open=[ParallelOpenFilePath]");
ParallelDataFile=getImageID;
BaseFileName=substring(File.nameWithoutExtension,0,indexOf(File.nameWithoutExtension,"_Parallel_raw_values"));
getDimensions(width, height, channels, slices, frames);
dir0=File.directory;

time1ArrayPar=newArray(height);
expt1ArrayPar=newArray(height);
time2ArrayPar=newArray(height);
expt2ArrayPar=newArray(height);
time3ArrayPar=newArray(height);
expt3ArrayPar=newArray(height);
time4ArrayPar=newArray(height);
expt4ArrayPar=newArray(height);
time5ArrayPar=newArray(height);
expt5ArrayPar=newArray(height);
time6ArrayPar=newArray(height);
expt6ArrayPar=newArray(height);
time7ArrayPar=newArray(height);
expt7ArrayPar=newArray(height);
time8ArrayPar=newArray(height);
expt8ArrayPar=newArray(height);
time9ArrayPar=newArray(height);
expt9ArrayPar=newArray(height);
time10ArrayPar=newArray(height);
expt10ArrayPar=newArray(height);



selectImage(ParallelDataFile);
for(i=0;i<height;i++){
if(width>=4){
time1ArrayPar[i]=getPixel(2,i);
expt1ArrayPar[i]=getPixel(3,i);
}
if(width>=6){
time2ArrayPar[i]=getPixel(4,i);
expt2ArrayPar[i]=getPixel(5,i);
}
if(width>=8){
time3ArrayPar[i]=getPixel(6,i);
expt3ArrayPar[i]=getPixel(7,i);
}
if(width>=10){
time4ArrayPar[i]=getPixel(8,i);
expt4ArrayPar[i]=getPixel(9,i);
}
if(width>=12){
time5ArrayPar[i]=getPixel(10,i);
expt5ArrayPar[i]=getPixel(11,i);
}
if(width>=14){
time6ArrayPar[i]=getPixel(12,i);
expt6ArrayPar[i]=getPixel(13,i);
}
if(width>=16){
time7ArrayPar[i]=getPixel(14,i);
expt7ArrayPar[i]=getPixel(15,i);
}
if(width>=18){
time8ArrayPar[i]=getPixel(16,i);
expt8ArrayPar[i]=getPixel(17,i);
}
if(width>=20){
time9ArrayPar[i]=getPixel(18,i);
expt9ArrayPar[i]=getPixel(19,i);
}
if(width>=22){
time10ArrayPar[i]=getPixel(20,i);
expt10ArrayPar[i]=getPixel(21,i);
}
}
close();
//*******end parallel data into arrays****



//******************** put the Perpendicular data into arrays 
run("Text Image... ", "open=[PerpendicularOpenFilePath]");
PerpendicularDataFile=getImageID;
//BaseFileName=File.nameWithoutExtension;
getDimensions(width, height, channels, slices, frames);

time1ArrayPer=newArray(height);
expt1ArrayPer=newArray(height);
time2ArrayPer=newArray(height);
expt2ArrayPer=newArray(height);
time3ArrayPer=newArray(height);
expt3ArrayPer=newArray(height);
time4ArrayPer=newArray(height);
expt4ArrayPer=newArray(height);
time5ArrayPer=newArray(height);
expt5ArrayPer=newArray(height);
time6ArrayPer=newArray(height);
expt6ArrayPer=newArray(height);
time7ArrayPer=newArray(height);
expt7ArrayPer=newArray(height);
time8ArrayPer=newArray(height);
expt8ArrayPer=newArray(height);
time9ArrayPer=newArray(height);
expt9ArrayPer=newArray(height);
time10ArrayPer=newArray(height);
expt10ArrayPer=newArray(height);



selectImage(PerpendicularDataFile);
for(i=0;i<height;i++){
if(width>=4){
time1ArrayPer[i]=getPixel(2,i);
expt1ArrayPer[i]=getPixel(3,i);
}
if(width>=6){
time2ArrayPer[i]=getPixel(4,i);
expt2ArrayPer[i]=getPixel(5,i);
}
if(width>=8){
time3ArrayPer[i]=getPixel(6,i);
expt3ArrayPer[i]=getPixel(7,i);
}
if(width>=10){
time4ArrayPer[i]=getPixel(8,i);
expt4ArrayPer[i]=getPixel(9,i);
}
if(width>=12){
time5ArrayPer[i]=getPixel(10,i);
expt5ArrayPer[i]=getPixel(11,i);
}
if(width>=14){
time6ArrayPer[i]=getPixel(12,i);
expt6ArrayPer[i]=getPixel(13,i);
}
if(width>=16){
time7ArrayPer[i]=getPixel(14,i);
expt7ArrayPer[i]=getPixel(15,i);
}
if(width>=18){
time8ArrayPer[i]=getPixel(16,i);
expt8ArrayPer[i]=getPixel(17,i);
}
if(width>=20){
time9ArrayPer[i]=getPixel(18,i);
expt9ArrayPer[i]=getPixel(19,i);
}
if(width>=22){
time10ArrayPer[i]=getPixel(20,i);
expt10ArrayPer[i]=getPixel(21,i);
}
}
close();
//*******end Perpendicular data into arrays****







tempArrayOfArraysPar=newArray("expt1ArrayPar","expt2ArrayPar","expt3ArrayPar","expt4ArrayPar","expt5ArrayPar","expt6ArrayPar","expt7ArrayPar","expt8ArrayPar","expt9ArrayPar","expt10ArrayPar");
tempArrayOfArraysPer=newArray("expt1ArrayPer","expt2ArrayPer","expt3ArrayPer","expt4ArrayPer","expt5ArrayPer","expt6ArrayPer","expt7ArrayPer","expt8ArrayPer","expt9ArrayPer","expt10ArrayPer");
tempArrayOfTimeArraysPar=newArray("time1ArrayPar","time2ArrayPar","time3ArrayPar","time4ArrayPar","time5ArrayPar","time6ArrayPar","time7ArrayPar","time8ArrayPar","time9ArrayPar","time10ArrayPar");
ArrayOfArraysPar=newArray(width/2-1);
ArrayOfArraysPer=newArray(width/2-1);
ArrayOfTimeArraysPar=newArray(width/2-1);
for(i=0;i<width/2-1;i++){
ArrayOfArraysPar[i]=tempArrayOfArraysPar[i];
ArrayOfArraysPer[i]=tempArrayOfArraysPer[i];
ArrayOfTimeArraysPar[i]=tempArrayOfTimeArraysPar[i];
}



decayFitParametersArray=newArray("expt"+","+"gFactor"+","+"CycleNumber"+","+"Parallel A"+","+"Parallel k"+","+"Parallel offset"+","+"Parallel R2"+","+"Perpendicular A"+","+"Perpendicular k"+","+"Perpendicular offset"+","+"Perpendicular R2"+","+"TotalF A"+","+"TotalF k"+","+"TotalF offset"+","+"TotalF R2");
RSFitParametersArray=newArray("expt"+","+"gFactor"+","+"CycleNumber"+","+"R0"+","+"NAObj"+","+"indexOfRefraction"+","+"ret0"+","+"slope/delta_r"+","+"ret1"+","+"R2"+","+"SE"+","+"ret0_Error"+","+"delta_r_Error"+","+"ret1_Error");

exptArrayR=newArray(height);//Anisotropy
exptArrayF=newArray(height);//total fluorescence

	Dialog.create("Input cycle information");
	Dialog.addMessage("*****************************************************************************");
	Dialog.addNumber("Number_of_cycles",3);
	Dialog.addNumber("Images_per_cycle",300);
	Dialog.addMessage("");
	Dialog.show();
	
numOfCycles=Dialog.getNumber();
numImagesPerCycle=Dialog.getNumber();
	


for(i=0;i<ArrayOfArraysPar.length;i++){
exptNum=IJ.pad(i+1,2);

if(i==0){
	extractAnisotropyData(dir0, exptNum, expt1ArrayPar, expt1ArrayPer, gFactor, NAObj, indexRefraction, time1ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt1ArrayPar, expt1ArrayPer, gFactor, time1ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt1ArrayPar, expt1ArrayPer, gFactor, R0, NAObj, indexRefraction, time1ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==1){
	extractAnisotropyData(dir0, exptNum, expt2ArrayPar, expt2ArrayPer, gFactor, NAObj, indexRefraction, time2ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt2ArrayPar, expt2ArrayPer, gFactor, time2ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt2ArrayPar, expt2ArrayPer, gFactor, R0, NAObj, indexRefraction, time2ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==2){
	extractAnisotropyData(dir0, exptNum, expt3ArrayPar, expt3ArrayPer, gFactor, NAObj, indexRefraction, time3ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt3ArrayPar, expt3ArrayPer, gFactor, time3ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt3ArrayPar, expt3ArrayPer, gFactor, R0, NAObj, indexRefraction, time3ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==3){
	extractAnisotropyData(dir0, exptNum, expt4ArrayPar, expt4ArrayPer, gFactor, NAObj, indexRefraction, time4ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt4ArrayPar, expt4ArrayPer, gFactor, time4ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt4ArrayPar, expt4ArrayPer, gFactor, R0, NAObj, indexRefraction, time4ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==4){
	extractAnisotropyData(dir0, exptNum, expt5ArrayPar, expt5ArrayPer, gFactor, NAObj, indexRefraction, time5ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt5ArrayPar, expt5ArrayPer, gFactor, time5ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt5ArrayPar, expt5ArrayPer, gFactor, R0, NAObj, indexRefraction, time5ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==5){
	extractAnisotropyData(dir0, exptNum, expt6ArrayPar, expt6ArrayPer, gFactor, NAObj, indexRefraction, time6ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt6ArrayPar, expt6ArrayPer, gFactor, time6ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt6ArrayPar, expt6ArrayPer, gFactor, R0, NAObj, indexRefraction, time6ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==6){
	extractAnisotropyData(dir0, exptNum, expt7ArrayPar, expt7ArrayPer, gFactor, NAObj, indexRefraction, time7ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt7ArrayPar, expt7ArrayPer, gFactor, time7ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt7ArrayPar, expt7ArrayPer, gFactor, R0, NAObj, indexRefraction, time7ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==7){
	extractAnisotropyData(dir0, exptNum, expt8ArrayPar, expt8ArrayPer, gFactor, NAObj, indexRefraction, time8ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt8ArrayPar, expt8ArrayPer, gFactor, time8ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt8ArrayPar, expt8ArrayPer, gFactor, R0, NAObj, indexRefraction, time8ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==8){
	extractAnisotropyData(dir0, exptNum, expt9ArrayPar, expt9ArrayPer, gFactor, NAObj, indexRefraction, time9ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt9ArrayPar, expt9ArrayPer, gFactor, time9ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt9ArrayPar, expt9ArrayPer, gFactor, R0, NAObj, indexRefraction, time9ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}

if(i==9){
	extractAnisotropyData(dir0, exptNum, expt10ArrayPar, expt10ArrayPer, gFactor, NAObj, indexRefraction, time10ArrayPar);
	decayFitParametersArray=fitDecayAnisotropyData(dir0, exptNum, expt10ArrayPar, expt10ArrayPer, gFactor, time10ArrayPar, numOfCycles, numImagesPerCycle, decayFitParametersArray);
	RSFitParametersArray=fitAnisotropyData(dir0, exptNum, expt10ArrayPar, expt10ArrayPer, gFactor, R0, NAObj, indexRefraction, time10ArrayPar, fitFractionCutOff, numOfCycles, numImagesPerCycle, RSFitParametersArray);
}


}//for(i=0;i<ArrayOfArraysPar.length;i++)

decayFitParameters_data=File.open(dir0+BaseFileName+"_DecayFitParameters.txt");
for(k=0;k<decayFitParametersArray.length;k++){
print(decayFitParameters_data,decayFitParametersArray[k]);
}
File.close(decayFitParameters_data);

if(File.exists(dir0+BaseFileName+"_DecayFitParameters_data.csv")){
File.delete(dir0+BaseFileName+"_DecayFitParameters_data.csv");
File.rename(dir0+BaseFileName+"_DecayFitParameters.txt",dir0+BaseFileName+"_DecayFitParameters_data.csv");
}else{
File.rename(dir0+BaseFileName+"_DecayFitParameters.txt",dir0+BaseFileName+"_DecayFitParameters_data.csv");	
}

RSFitParameters_data=File.open(dir0+BaseFileName+"_RSFitParameters.txt");
for(k=0;k<RSFitParametersArray.length;k++){
print(RSFitParameters_data,RSFitParametersArray[k]);
}
File.close(RSFitParameters_data);

if(File.exists(dir0+BaseFileName+"_RSFitParameters.csv")){
File.delete(dir0+BaseFileName+"_RSFitParameters_data.csv");
File.rename(dir0+BaseFileName+"_RSFitParameters.txt",dir0+BaseFileName+"_RSFitParameters_data.csv");
}else{
File.rename(dir0+BaseFileName+"_RSFitParameters.txt",dir0+BaseFileName+"_RSFitParameters_data.csv");	
}

function fitAnisotropyData(dir, exptN, exptArrayPar, exptArrayPer, gFac, R0, ObjNA, refractiveIndex, timeArray, fractionCutOffToFit, numCycles, imagesPerCycle, arrayToReturn){

for(i=0;i<numCycles;i++){
	x=Array.slice(timeArray,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	Array.getStatistics(x,xminVal,xmaxVal,xmeanVal,xsdVal);
	x=arraySubtractValue(x,xminVal);
	yPar=Array.slice(exptArrayPar,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	yPer=Array.slice(exptArrayPer,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	exptArray2GPer=Array.copy(yPer);
	exptArray2GPer=arrayMultiplyByValue(exptArray2GPer,2*gFac);
	exptArrayF=addTwoArrays(yPar,exptArray2GPer);
	Array.getStatistics(exptArrayF, min, max, mean, stdDev);
	exptArrayN=Array.copy(exptArrayF);
	exptArrayN=arrayDivideByValue(exptArrayN,max);
	exptArrayNPS=Array.copy(exptArrayN);
	exptArrayNPS=valueSubtractArray(1,exptArrayNPS);
	exptArrayR=calculateAnisotropyTwoArrays(yPar,yPer,gFac);

	if(calculateUncorrected){
		exptArrayCR=calculateAnisotropyTwoArrays(yPar,yPer,gFac);
	}else{
		exptArrayCR=calculateCorrectedAnisotropyTwoArrays(yPar,yPer,gFac,ObjNA,refractiveIndex);
	}
	
	fitCutOff=findFitFractionCutoff(exptArrayNPS,fractionCutOffToFit);
	if(fitCutOff==0){
		ret=NaN;
		delta_r=NaN;//the slope
		r0=NaN;
		R2=NaN;
	}else{
		exptArrayNToFit=Array.trim(exptArrayNPS,fitCutOff);
		exptArrayRToFit=Array.trim(exptArrayCR,fitCutOff);
		
		//pre-fit to linear to get r1
		m=(exptArrayRToFit[exptArrayRToFit.length-1]-exptArrayRToFit[0])/(exptArrayNToFit[exptArrayNToFit.length-1]-exptArrayNToFit[0]);
		initialGuessesL=newArray(0.3,m);//a,b of y = a + bx;
		Fit.doFit("Straight Line",exptArrayNToFit,exptArrayRToFit,initialGuessesL);
		//Fit.plot;
		ret=Fit.p(0);
		delta_r=Fit.p(1);//the slope
		r0=ret+delta_r;
		R2=Fit.rSquared;
		sumOfX=0;
		sumOfSqrX=0;
		theFitArray=newArray(exptArrayNToFit.length);
		for(tf=0;tf<exptArrayNToFit.length;tf++){
			theFitArray[tf]=Fit.f(exptArrayNToFit[tf]);
			sumOfX=sumOfX+exptArrayNToFit[tf];
			sumOfSqrX=sumOfSqrX+(exptArrayNToFit[tf]*exptArrayNToFit[tf]);
		}
		SE=calculateSERegression(exptArrayRToFit,theFitArray);
		delta_r_Error=SE*(sqrt(exptArrayNToFit.length/((exptArrayNToFit.length*sumOfSqrX)-(sumOfX*sumOfX))));
		ret_Error=SE*(sqrt(sumOfSqrX/((exptArrayNToFit.length*sumOfSqrX)-(sumOfX*sumOfX))));
		r0_Error=sqrt((delta_r_Error*delta_r_Error)+(ret_Error*ret_Error));
	}

	fitParameters=toString(exptN+","+gFac+","+i+","+R0+","+ObjNA+","+refractiveIndex+","+r0+","+delta_r+","+ret+","+R2+","+SE+","+r0_Error+","+delta_r_Error+","+ret_Error);
	arrayToReturn=appendToArray(fitParameters,arrayToReturn);
}
return arrayToReturn;
}//end function fitAnisotropyData


function fitDecayAnisotropyData(dir, exptN, exptArrayPar, exptArrayPer, gFac, timeArray, numCycles, imagesPerCycle, arrayToReturn){

for(i=0;i<numCycles;i++){
	x=Array.slice(timeArray,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	Array.getStatistics(x,xminVal,xmaxVal,xmeanVal,xsdVal);
	x=arraySubtractValue(x,xminVal);
	yPar=Array.slice(exptArrayPar,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	yPer=Array.slice(exptArrayPer,i*imagesPerCycle,(i*imagesPerCycle)+(imagesPerCycle-1));
	exptArray2GPer=Array.copy(yPer);
	exptArray2GPer=arrayMultiplyByValue(exptArray2GPer,2*gFac);
	exptArrayF=addTwoArrays(yPar,exptArray2GPer);
	Array.getStatistics(exptArrayF, min, max, mean, stdDev);
	exptArrayN=Array.copy(exptArrayF);
	exptArrayN=arrayDivideByValue(exptArrayN,max);
	exptArrayR=calculateAnisotropyTwoArrays(yPar,yPer,gFac);
	
	single_exp_fit="y = a*exp(-b*x)+c";
	//Fit Parallel
	initialGuessesPar=newArray(yPar[0],1,yPar[yPar.length-1]);
	Fit.doFit(single_exp_fit,x,yPar,initialGuessesPar);
	//Fit.plot;
	aPar=Fit.p(0);
	bPar=Fit.p(1);
	cPar=Fit.p(2);
	RPar=Fit.rSquared;

	//Fit Perpendicular
	initialGuessesPer=newArray(yPer[0],1,yPer[yPer.length-1]);
	Fit.doFit(single_exp_fit,x,yPer,initialGuessesPer);
	//Fit.plot;
	aPer=Fit.p(0);
	bPer=Fit.p(1);
	cPer=Fit.p(2);
	RPer=Fit.rSquared;

	//Fit Total Fluorescence
	initialGuessesF=newArray(exptArrayF[0],1,exptArrayF[exptArrayF.length-1]);
	Fit.doFit(single_exp_fit,x,exptArrayF,initialGuessesF);
	//Fit.plot;
	aF=Fit.p(0);
	bF=Fit.p(1);
	cF=Fit.p(2);
	RF=Fit.rSquared;
	fitParameters=exptN+","+gFac+","+i+","+aPar+","+bPar+","+cPar+","+RPar+","+aPer+","+bPer+","+cPer+","+RPer+","+aF+","+bF+","+cF+","+RF;
	arrayToReturn=appendToArray(fitParameters,arrayToReturn);
}
return arrayToReturn;
}//end function fitAnisotropyData






function findFitFractionCutoff(theArray,theCutOff){
cutOffToReturn=0;
for(i=0;i<theArray.length-1;i++){
	if(theArray[i+1]>theCutOff && theArray[i]<=theCutOff)
	cutOffToReturn=i;
}
return cutOffToReturn;
}

function extractAnisotropyData(dir, exptN, exptArrayPar, exptArrayPer, gFac, ObjNA, refractiveIndex, timeArray){
	anisotropy_data=File.open(dir+"expt"+exptN+"_anisotropy.txt");
	print(anisotropy_data, "expt"+exptN+"_Time,"+"expt"+exptN+"_gFactor,"+"expt"+exptN+"_ParallelF,"+"expt"+exptN+"_PerpendicularF,"+"expt"+exptN+"_TotalF,"+"expt"+exptN+"_normalF,"+"expt"+exptN+"_Fluorophore_Photoswitched,"+"expt"+exptN+"_r,"+"expt"+exptN+"_corrected_r");
	exptArray2GPer=Array.copy(exptArrayPer);
	exptArray2GPer=arrayMultiplyByValue(exptArray2GPer,2*gFac);
	exptArrayF=addTwoArrays(exptArrayPar,exptArray2GPer);
	Array.getStatistics(exptArrayF, min, max, mean, stdDev);
	exptArrayN=Array.copy(exptArrayF);
	exptArrayN=arrayDivideByValue(exptArrayN,max);
	exptArrayNPS=Array.copy(exptArrayN);
	exptArrayNPS=valueSubtractArray(1,exptArrayNPS);
	exptArrayR=calculateAnisotropyTwoArrays(exptArrayPar,exptArrayPer,gFac);
	exptArrayCR=calculateCorrectedAnisotropyTwoArrays(exptArrayPar,exptArrayPer,gFac,ObjNA,refractiveIndex);
	for(j=0;j<exptArrayPar.length;j++){
	print(anisotropy_data, timeArray[j]+","+gFac+","+exptArrayPar[j]+","+exptArrayPer[j]+","+exptArrayF[j]+","+exptArrayN[j]+","+exptArrayNPS[j]+","+exptArrayR[j]+","+exptArrayCR[j]);
	}
	File.close(anisotropy_data);

	if(File.exists(dir+"expt"+exptN+"_anisotropy_data.csv")){
	File.delete(dir+"expt"+exptN+"_anisotropy_data.csv");
	File.rename(dir+"expt"+exptN+"_anisotropy.txt",dir+"expt"+exptN+"_anisotropy_data.csv");
	}else{
	File.rename(dir+"expt"+exptN+"_anisotropy.txt",dir+"expt"+exptN+"_anisotropy_data.csv");	
	}
}

function valueSubtractArray(value,theArray){
for(i=0;i<theArray.length;i++){
theArray[i]=value-theArray[i];
}
return theArray;
}

function arraySubtractValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]-value;
}
return theArray;
}


function arrayAddValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]+value;
}
return theArray;
}

function arrayDivideByValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]/value;
}
return theArray;
}

function arrayMultiplyByValue(theArray,value){
for(i=0;i<theArray.length;i++){
theArray[i]=theArray[i]*value;
}
return theArray;
}

function addTwoArrays(A1,A2){
	if(A1.length==A2.length){
		sumArray=newArray(A1.length);
		for(i=0;i<A1.length;i++){
			sumArray[i]=A1[i]+A2[i];
		}
	}
	return sumArray;
}

function subtractTwoArrays(A1,A2){
	if(A1.length==A2.length){
		sumArray=newArray(A1.length);
		for(i=0;i<A1.length;i++){
			sumArray[i]=A1[i]-A2[i];
		}
	}
	return sumArray;
}


function multiplyTwoArrays(A1,A2){
	if(A1.length==A2.length){
		sumArray=newArray(A1.length);
		for(i=0;i<A1.length;i++){
			sumArray[i]=A1[i]*A2[i];
		}
	}
	return sumArray;
}


function divideTwoArrays(A1,A2){
	if(A1.length==A2.length){
		sumArray=newArray(A1.length);
		for(i=0;i<A1.length;i++){
			sumArray[i]=A1[i]/A2[i];
		}
	}
	return sumArray;
}


function calculateAnisotropyTwoArrays(ParallelArray,PerpendicularArray,g_Factor){
	if(ParallelArray.length==PerpendicularArray.length){
		sumArray=newArray(ParallelArray.length);
		for(i=0;i<ParallelArray.length;i++){
			sumArray[i]=(ParallelArray[i]-(g_Factor*PerpendicularArray[i]))/(ParallelArray[i]+(2*g_Factor*PerpendicularArray[i]));
		}
	}
	return sumArray;
}

function calculateCorrectedAnisotropyTwoArrays(ParallelArray,PerpendicularArray,g_Factor,NA,n){
	alpha=asin(NA/n);
	//print("alpha="+alpha);
	Ka=1/3*(2-(3*cos(alpha))+(pow(cos(alpha),3)));
	Kb=1/12*(1-(3*cos(alpha))+(3*(pow(cos(alpha),2)))-(pow(cos(alpha),3)));
	Kc=1/4*(5-(3*cos(alpha))-(pow(cos(alpha),2))-(pow(cos(alpha),3)));
	//print("Ka="+Ka+"   Kb="+Kb+"   Kc="+Kc);
	
	if(ParallelArray.length==PerpendicularArray.length){
		sumArray=newArray(ParallelArray.length);
		for(i=0;i<ParallelArray.length;i++){
			Iy=((ParallelArray[i]*Kb)-(g_Factor*PerpendicularArray[i]*Kc))/((Kb*Kb)+(Ka*Kb)-(Kc*Kc)-(Ka*Kc));
			Iz=((g_Factor*PerpendicularArray[i]*Kb)+(g_Factor*PerpendicularArray[i]*Ka)-(ParallelArray[i]*Kc)-(ParallelArray[i]*Ka))/((Kb*Kb)+(Ka*Kb)-(Kc*Kc)-(Ka*Kc));		
			sumArray[i]=(Iz-Iy)/(Iz+(2*Iy));
		}
	}
	return sumArray;
}

function appendToArray(value, array) {
temparray=newArray(lengthOf(array)+1);
for(i=0;i<lengthOf(array);i++){
temparray[i]=array[i];
}
temparray[lengthOf(temparray)-1]=value;
array=temparray;
return array;
}

function findStringFragment(theString,theFragment){
	stringFound=false;
for(i=0;i<=lengthOf(theString)-lengthOf(theFragment);i++){
	check4Frag=substring(theString,i,i+lengthOf(theFragment));
	if(check4Frag==theFragment)
	stringFound=true;
}
	return stringFound;
}


function  calculateSERegression(dataArray,fitArray){
	theSum=0;
	for(i=0;i<dataArray.length;i++){
		theSum=theSum+pow((fitArray[i]-dataArray[i]),2);
	}
	StdErr=sqrt(theSum/(dataArray.length-2-1));
	return StdErr;
}

}//end of macro