#pragma rtGlobals=1		// Use modern global access method.

Function ZAPSimAnal(InputWave)
	wave InputWave
	 
	String OutputWaveName = NameOfWave(InputWave) + "Z" 
	String SmoothWaveName = NameOfWave(InputWave) + "ZS" 
	
	Make /O/N=625 $outputwavename
	Wave outputwave=$outputwavename
	Variable ZAPAmp = 10
	make /o /N=1000000 ZapI
	SetScale/I x 0,25,"", ZapI
	ZapI=ZAPAmp*sin(2*3.141592654*((25*x)/(2*25))*x)
	FFT /dest=ZapIFFT zapi
	FFT /dest=zapvfft /RP=[2000,1001999] InputWave
	make /o /N=1024 /C ZapImp
	ZapImp=ZapVFFT/ZapIFFT
	OutputWave=sqrt(magsqr(zapimp))
	OutputWave[0]=Outputwave[3]
	OutputWave[1]=Outputwave[3]
	OutputWave[2]=Outputwave[3]
	OutputWave*=1000
	SetScale/I x 0,25,"", OutputWave
	Duplicate /O outputwave, $smoothwavename
	smooth 25, $smoothwavename
	wavestats /q /R=(0,25) $smoothwavename
	killwaves zapvfft, zapifft
	Duplicate /O  $smoothwavename, outputwave
	
	duplicate /o Inputwave, ZAPVoltage
	deletepoints 0,50, ZAPVoltage
	SetScale/I x 0,25,"", ZapVoltage

	print "fR =", v_maxloc
	print "Q =", outputwave(V_maxloc)/outputwave(0.5)
	print "Zmax =", v_max
	make /o /n=(3,1) returnwave
	returnwave[0]=v_maxloc
	returnwave[1]=outputwave(V_maxloc)/outputwave(0.5)
	returnwave[2]=v_max
	return returnwave
End	

Function ZAPSimFamilyAnal(InputBaseName, outputwavename)
	String InputBaseName
	String Outputwavename
	
	String wavename1
	String fRwavename
	String Qwavename
	String Zmaxwavename
	
	String Filename

	Variable i
	Variable count=0
	
	sprintf fRwavename, "%s_fR", Outputwavename
	sprintf Qwavename, "%s_Q", Outputwavename
	sprintf Zmaxwavename, "%s_Zmax", Outputwavename
		
	Make /O/N=500 $fRwavename
	Make /O/N=500 $Qwavename
	Make /O/N=500 $Zmaxwavename

	make /o/n=(3,1) rwave
	
	Wave fRwave=$fRwavename
	Wave QWave=$Qwavename
	Wave ZmaxWave=$Zmaxwavename
	
	//for(i=0.1; i<=3.0; i+=0.1)
	//for(i=80; i>=55; i-=5)	
	for(i=1; i<=68; i+=1)	
		
	//	sprintf Filename, "Macintosh HD:Users:manisha:FUN:LFP_Results:Osc8:Transfer:Zap_gh_0.000028_apical[65]_ZAPAmp_%1.6f.txt",i
	//sprintf Filename, "Macintosh HD:Users:manisha:FUN:LFP_Results:HopefullyFinalGrad:CNL_ParameterValidation:Fig_Zap_%s_gh160_%d.txt",inputbasename,i
	//sprintf Filename, "Macintosh HD:Users:manisha:FUN:LFP_Results:newGradient:Osc8:Trial:Ra_50_15_Rm_55k_10k:Zap_%s_%d.txt",inputbasename,i
	//sprintf Filename, "Macintosh HD:Users:manisha:FUN!!!!:LFP model:OscFreq5:MorphologicallyRealistic:ParameterValidation:Zap_%s_%d.txt",inputbasename,i
	//for(i=1/10; i<=2; i+=1/100)
		//sprintf Filename, "thodiHD:Users:rishi:Work:Influence Field:Output:%s%4.2f.txt",InputBaseName,i
		//sprintf Filename, "thodiHD:Users:rishi:Work:Influence Field:Output:%s%d.txt",InputBaseName,i
		//sprintf Filename, "RishiHD:Users:rishi:Work:Depletion:Simulations:Output:%s%d.txt",InputBaseName,i
		
		
		//sprintf Filename, "Macintosh HD:Users:apple:FUN:PNAS_newResults:NewMorph_reloaded:ZapFinal:Zap_%d.txt",i
		sprintf Filename, "Macintosh HD:Users:apple:FUN:PNAS_newResults:NewMorph_reloaded:ZapFinal:Zap_Trans_%d.txt",i
		
		print i, filename
		loadwave /J/N filename		
		sprintf wavename1, "Zap_%s_%d",outputwavename,i
		duplicate /O wave0, $wavename1

		print "THIS:" , wavename1
		rwave=ZAPSimAnal($wavename1)
		ZAPPhase($wavename1)
		
		fRwave[count]=rwave[0]
		Qwave[count]=rwave[1]
		Zmaxwave[count]=rwave[2]
		print rwave
		count=count+1
	endfor
	redimension /N=(count) fRwave
	redimension /N=(count) Qwave
	redimension /N=(count) Zmaxwave
End	

Function ZAPPhase(InputWave)
	wave InputWave
	 
	String OutputWaveName = NameOfWave(InputWave) + "P" 
	String SmoothWaveName = NameOfWave(InputWave) + "PS" 
	
	Make /O/N=625 $outputwavename
	Wave outputwave=$outputwavename

	Variable ZAPAmp = 10
	make /o /N=1000000 ZapI
	SetScale/I x 0,25,"", ZapI
	ZapI=ZAPAmp*sin(2*3.141592654*((25*x)/(2*25))*x)
	FFT /dest=ZapcurrentFFT zapi
	FFT /dest=zapvoltagefft /RP=[2000,1001999] InputWave
	make /o /N=1024 /C ZapImpedance
	ZapImpedance=ZapVoltageFFT/ZapcurrentFFT
	OutputWave=-imag(r2polar(zapimpedance))*180/pi
	OutputWave[0]=Outputwave[3]
	OutputWave[1]=Outputwave[3]
	OutputWave[2]=Outputwave[3]
 	SetScale/I x 0,25,"", OutputWave
 	Duplicate /O outputwave, $smoothwavename
	smooth 25, $smoothwavename
End	

// Designed for a ZAP20/ZAP15 stimulus.

Function ZAPExperimentPhase(InputWave)
	wave InputWave
	 
	String OutputWaveName = NameOfWave(InputWave) + "P" 
	String SmoothWaveName = NameOfWave(InputWave) + "PS" 

	Variable ZAPDur
	Variable ZAPMaxFreq
	Variable InputDim
	
	InputDim=numpnts(InputWave) 
	
	if (InputDim==215000)
		ZAPDur=20000
		ZAPMaxFreq=20
	elseif (InputDim==165000)
		ZAPDur=15000
		ZAPMaxFreq=15		
	endif
	
	make /o/n=(1) ZapCurrent					// The ZAP input wave
	make /o/n=(1) /C ZAPImp					// Complex impedence wave
	duplicate /o InputWave ZapVoltage		// The ZAP input wave
	deletepoints 0,15000,ZAPVoltage
		
	Make /O/N=(1) $outputwavename			// Output wave.
	Wave ZapOutput=$outputwavename
	Wave ZapSmooth=$smoothwavename

	Redimension /n=512  zapimp
	Redimension/N=(ZAPMaxFreq*ZAPMaxFreq) ZAPOutput
	Redimension/N=(InputDim-15000) ZAPCurrent
	SetScale/I x 0,ZAPDur,"", ZapCurrent, ZAPVoltage


	ZapCurrent=50*sin(2*3.141592654*((ZAPMaxFreq*x)/(2*ZAPDur))*x/1000)
	
	FFT /dest=ZapcurrentFFT ZapCurrent
	FFT /dest=zapvoltagefft ZAPVoltage

	ZapImp=ZapVoltageFFT/ZapcurrentFFT
	ZapOutput=-imag(r2polar(ZapImp))*180/pi
	SetScale/I x 0, ZAPMaxFreq,"", ZapOutput
	Duplicate /O zapoutput, $smoothwavename
	smooth 25, $smoothwavename
	ZapOutput[0,5]=nan
	zapsmooth[0,5]=nan
End	

Function InductiveArea (InputWave)
	wave InputWave
	
	Variable i
	String PhaseWaveName = NameOfWave(InputWave) + "P" 
	print PhaseWaveName

	ZAPPhase(Inputwave)
	Wave PhaseWave=$Phasewavename
	
	duplicate /o PhaseWave, AnalysisWave
	
	AnalysisWave= (PhaseWave[p]<0)?0:Phasewave[p]
	
	return area (Analysiswave)
End	
	
Function InductiveAreaFamilyAnal (InputBaseName, outputwavename)
	String InputBaseName
	String Outputwavename
	
	String AUCwavename, PSWavename
	
	String Filename

	Variable i
	Variable count=0
	
	sprintf AUCwavename, "%s_AUC", Outputwavename
		
	Make /O/N=500 $AUCwavename
	 	
	Wave AUCWave=$AUCwavename
	
	
	//for(i=-130; i<=-50; i+=1)	
	//for(i=0.1; i<=3; i+=0.1)
	//for(i=1; i<=25; i+=1)	
	for(i=0.05; i<=0.3; i+=0.05)	
		sprintf Filename, "Macintosh HD:Users:manisha:FUN:LFP_Results:Osc8:Transfer:Zap_gh_0.000028_apical[65]_ZAPAmp_%1.6f.txt",i
	//for(i=1/10; i<=2; i+=1/100)
		//sprintf Filename, "RishiHD:Users:rishi:Oscillations:Simulations:Single:Output:%s%4.3f.txt",InputBaseName,i
		//sprintf Filename, "RishiHD:Users:rishi:Oscillations:Simulations:Multi:Output:%s%d.txt",InputBaseName,i
		//sprintf Filename, "RishiHD:Users:rishi:Oscillations:Simulations:Single:Output:%s%d.txt",InputBaseName,i
		print i, filename
		loadwave /J/N filename
		sprintf PSWavename, "%s%d",outputwavename, i*1000
		duplicate /O wave0, $PSWavename
		Wave PSWave=$PSWavename
		AUCwave[count]=InductiveArea(PSWave)
		print AUCwave[count]
		count=count+1
	endfor
	redimension /N=(count) AUCwave
	AUCwave*=pi/180 // Convert to Radians.
End	

//	for (i=numpnts(Inputwave)-1; i>=0; i=i-1)
//		if (AnalysisWave[i]!=0)
//			break;
//		endif
//	endfor
//	
//	DeletePoints i+2, numpnts(Inputwave)-i-1, Analysiswave

Function ZAPSimAnal35(InputWave)
	wave InputWave
	 
	String OutputWaveName = NameOfWave(InputWave) + "Z" 
	String SmoothWaveName = NameOfWave(InputWave) + "ZS" 
	
	Make /O/N=1225 $outputwavename
	Wave outputwave=$outputwavename
	
	make /o /N=1400000 ZapI
	SetScale/I x 0,35,"", ZapI
	ZapI=50*sin(2*3.141592654*((35*x)/(2*35))*x)
	FFT /dest=ZapIFFT zapi
	FFT /dest=zapvfft /RP=[2000,1401999] InputWave
	make /o /N=2048 /C ZapImp
	ZapImp=ZapVFFT/ZapIFFT
	OutputWave=sqrt(magsqr(zapimp))
	OutputWave[0]=Outputwave[3]
	OutputWave[1]=Outputwave[3]
	OutputWave[2]=Outputwave[3]
	OutputWave*=1000
	SetScale/I x 0,35,"", OutputWave
	Duplicate /O outputwave, $smoothwavename
	smooth 25, $smoothwavename
	wavestats /q /R=(0,35) $smoothwavename
	killwaves zapvfft, zapifft
	Duplicate /O  $smoothwavename, outputwave
	
	duplicate /o Inputwave, ZAPVoltage
	deletepoints 0,50, ZAPVoltage
	SetScale/I x 0,35,"", ZapVoltage

	print "fR =", v_maxloc
	print "Q =" outputwave(V_maxloc)/outputwave(0.5)
	print "Zmax =", v_max
	make /o /n=(3,1) returnwave
	returnwave[0]=v_maxloc
	returnwave[1]=outputwave(V_maxloc)/outputwave(0.5)
	returnwave[2]=v_max
	return returnwave
End