#pragma rtGlobals=1		// Use modern global access method.
#include <differentiate XY>
#include <Median>



Window Panel01() : Panel  
	PauseUpdate; Silent 1		// building window...
	VARIABLE/G NPFSTime, NPFEtime
	NewPanel /W=(356,388,1025,883)
	ShowTools
	SetDrawLayer UserBack
	SetDrawEnv linethick= 0
	DrawRect 26,12,563,48
	SetDrawEnv fname= "Century Gothic",fsize= 25,fstyle= 1 
	DrawText 32,43,"PARGAN (PARticle Growth And Nucleation)"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif"
	DrawText 188,101,"Time wave"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif"
	DrawText 187,176,"dN/dLogDp (cm-3)"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif",fstyle= 1
	DrawText 30,372,"Calculate Growth Rate (nm.h-1)"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif",fstyle= 1
	DrawText 294,373,"Calculate Nucleation Rate (cm-3.s-1)"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif",fstyle= 1
	DrawText 28,74,"Input files "
	SetDrawEnv linethick= 3,fname= "MS Sans Serif"
	DrawText 186,207,"Raw count matrix"
	SetDrawEnv linethick= 3,fname= "MS Sans Serif"
	DrawText 186,137,"Diameter wave (nm)"
	SetDrawEnv linefgc= (65280,0,0),fillfgc= (65280,0,0)
	DrawRect 22,77,184,223
	Button Obtaining_growth_rate,pos={27,377},size={167,30},proc=Growth_Rate,title="Calculate Growth Rate"
	Button Graph_growthvstime,pos={28,413},size={167,28},proc=Call_Graph_growth_vs_time,title="Plot:  Growth Rate vs Time"
	Button Obtaining_nucl_rate,pos={294,378},size={185,28},proc=Nucl_Rate,title="Calculate Nucleation Rate (J)"
	Button Graph_nuclrate,pos={294,446},size={186,28},proc=Nucl_Rate_Graphs,title="Plot:  Nucl Rate vs Time"
	Button nuclrate_median,pos={294,412},size={185,28},proc=Get_Nucl_Rate_median,title="Get 1 min. interval median J"
	Button load_n1,pos={26,229},size={155,30},proc=Makebanana_func,title="Make Banana Plot"
	Button load_n1,fStyle=0
	Button load_dndlogdpmatrix,pos={27,155},size={155,30},proc=Input_dNdlogDp_matrix,title="Load dN_dlogD"
	Button load_dndlogdpmatrix,fStyle=0
	Button load_rawcounts_matrix,pos={26,190},size={155,30},proc=input_counts,title="Load Raw Counts matrix"
	Button load_rawcounts_matrix,fStyle=0
	SetVariable NPFStartTime,pos={29,266},size={155,16},title="NPF starting point"
	SetVariable NPFStartTime,value= NPFSTime
	SetVariable NPFendingTime,pos={31,286},size={155,16},title="NPF ending point"
	SetVariable NPFendingTime,value= NPFEtime
	Button load_n2,pos={29,308},size={155,30},proc=Select_NPFtime,title="Select NPF time"
	Button load_n2,fStyle=0
	Button load_dndlogdpmatrix1,pos={26,80},size={155,30},proc=input_time,title="Load time"
	Button load_dndlogdpmatrix1,fStyle=0
	Button load_dndlogdpmatrix2,pos={26,118},size={155,30},proc=input_diam,title="Load Dp"
	Button load_dndlogdpmatrix2,fStyle=0
	Button Graph_nuclrate1,pos={493,446},size={130,28},proc=save_GR_Jnuc,title="Save GR and Jnuc"
	Button load_dndlogdpmatrix3,pos={302,157},size={155,30},proc=Input_dNdlogDp_matrix,title="Load dN_dlogD"
	Button load_dndlogdpmatrix3,fStyle=0
	Button load_rawcounts_matrix1,pos={303,190},size={155,30},proc=Input_cnt_matrix,title="Load Raw Counts matrix"
	Button load_rawcounts_matrix1,fStyle=0
EndMacro       
  
//SetVariable var_str_row_5,pos={309,382},size={174,20},title="Time Interval (min)"
//SetVariable var_str_row_5,help={"Enter the starting date"},value= K16
//Button nuclrate_median1,pos={489,377},size={102,28},proc=Get_Jmedian,title="Get median J"
 
NVAR NPFSTime, NPFEtime
    
//VARIABLE/G NPFSTime, NPFEtime
    
//FUNCTION Load_diffusionfactor(ctrlName) : ButtonControl    
//	STRING ctrlName 
//	STRING pathname=""  
//	STRING filename=""	
//	LOADWAVE/A/J/D/A=DC_factor/P=$pathname/K=0 filename
//	DUPLICATE/O DC_factor0 DC_factor  
//	KILLWAVES/Z DC_factor0
//END     

//FUNCTION input_time(ctrlName) : ButtonControl    
//	STRING ctrlName 
//	STRING pathname="" 
//	STRING filename=""	
//	LOADWAVE/A/J/D/A=t_wave/P=$pathname/K=0 filename
//	DUPLICATE/O t_wave0 t_wave  
//	KILLWAVES/Z t_wave0
//END    
//   
//FUNCTION input_diam(ctrlName) : ButtonControl                
//	STRING ctrlName 
//	STRING pathname=""
//	STRING filename=""	
//	LoadWave/A/J/D/A=diam/P=$pathname/K=0 filename 
//	DUPLICATE/O diam0 diam
//	MAKE/D/O/N=(NUMPNTS(diam)) radius
//	radius = diam/2. 
//	KILLWAVES/Z diam0
//END  

//FUNCTION Input_dNdlogDp_matrix(ctrlName) : ButtonControl    
//	STRING ctrlName 
//	STRING pathname=""
//	STRING filename=""	
//	LoadWave/J/M/U={0,0,0,0}/D/O/N=n/K=0/P=$pathname filename
//	DUPLICATE/O n0 n
//	KILLWAVES/Z n0  
//	
//	//############# Fitting lognormal fucntion to the data ####################
//	WAVE t_wave, diam
//	VARIABLE rows,cols
//	rows=DIMSIZE(t_wave,0)
//	cols=DIMSIZE(diam,0)  
//	MAKE/D/N=(rows,cols) dNdlogDp
//	dNdlogDp=n
//	//SMOOTH/B=3 2, dNdlogDp   
//	SMOOTH/B=2 5,dNdlogDp  
//	
//	MAKE/O/D/N=(3350) actual_N, smoothed_N
//	 VARIABLE i,j,nn,mm
//	 actual_N=NAN
//	 smoothed_N=NAN
//	 nn=0
//	 mm=0 
//	 FOR(i=0;i<rows;i+=1) 
//	  	FOR(j=0;j<cols;j+=1) 
//	  		actual_N[nn] = n[i][j]
//	  		smoothed_N[mm]=dNdlogDp[i][j]
//	  		nn=nn+1
//	  		mm=mm+1	  	
//	  	ENDFOR  
//	 ENDFOR  
////	VARIABLE i,j  
////	MAKE/O/D/N=(cols) dNdlogDp_wave,dNdlogDp_fitted_wave
////	MAKE/O/D/N=(rows,cols) dNdlogDp_fitted_matrix
////	MAKE/O/D/N=(rows) slope_lreg
////	MAKE/O/D/N=(2) W_coef
////	W_coef=NAN
////	slope_lreg=NAN
////	FOR(i=0;i<rows;i+=1) 
////      		FOR(j=0;j<cols;j+=1) 
////      			dNdlogDp_wave[j] = dNdlogDp[i][j]	        
////          	ENDFOR   
////		DUPLICATE/O/D dNdlogDp_wave dNdlogDp_wave_dum    
////		
////		CurveFit LogNormal  dNdlogDp_wave_dum /X=diam  /D=dNdlogDp_fitted_wave
////		CurveFit/M=2/W=0 line, dNdlogDp_fitted_wave/X=dNdlogDp_wave/D
////		slope_lreg[i]=W_coef[1]
////	
////		//CurveFit/M=2/W=0 poly 20,  dNdlogDp_wave_dum /X=diam  /D=dNdlogDp_fitted_wave
////		//CurveFit/NTHR=0 poly 7, dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
////		// CurveFit/NTHR=0 exp  dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
////		  
//////		display dNdlogDp_wave_dum vs diam
//////		AppendToGraph dNdlogDp_fitted_wave vs diam 
//////		ModifyGraph rgb(dNdlogDp_fitted_wave)=(63232,0,11520)
//////		ModIFyGraph log(bottom)=1 
//////		SetAxis bottom 3.2,1000	   
//////		ModifyGraph mode=4,rgb=(0,0,0) 	        
//////		PauseUpdate     
//////		String graphName   
//////		variable id=0
//////		graphName = WinName(id, 1)      
//////		if (strlen(graphName) == 0)
//////			break  
//////		endif    
//////		String fileNames
//////		SaveGraphCopy/P=$pathName/W=$graphName as fileNames
//////			  
////		FOR(j=0;j<cols;j+=1)         
////      			IF( dNdlogDp_fitted_wave[j] >=0.0)	 
////      			       dNdlogDp_fitted_matrix[i][j] = dNdlogDp_fitted_wave[j]	   
////      			ELSE
////      				dNdlogDp_fitted_matrix[i][j] =NAN
////      			ENDIF  
////          	ENDFOR   
////          	    
////	ENDFOR       
////	 
////	DUPLICATE/O/D dNdlogDp_fitted_matrix dNdlogDp              
//	
//	//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//	// Note that for dNdlogDp and raw counts matrix should be rows=rows_diam AND cols=rows_t_wave
//	//Since we have opposite formatting as required, we need to transpose dNdlogDp and raw count matrix's
//	 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//      MAKE/O/D/N=( (DIMSIZE(dNdlogDp,1)), (DIMSIZE(dNdlogDp,0)) ) dNdlogDp_dum
//      DUPLICATE/O/D dNdlogDp dNdlogDp_dum,dNdlogDp_plot
//      Redimension dNdlogDp_dum
//	MatrixTranspose dNdlogDp_dum
//	DUPLICATE/O/D dNdlogDp_dum dNdlogDp   
//	KILLWAVES dNdlogDp_dum   
//END 
//  

//// this input count matrix , same as dndlogdp format
//FUNCTION input_counts(ctrlName) : ButtonControl   
//	STRING ctrlName 
//	STRING pathname=""
//	STRING filename=""	
//	LoadWave/J/M/U={0,0,0,0}/D/O/N=cnt/K=0/P=$pathname filename
//	DUPLICATE/O cnt0 cnt
//	KILLWAVES/Z cnt0
//	
////	WAVE t_wave, diam
////	VARIABLE rows,cols,i,j
////	rows=DIMSIZE(t_wave,0)
////	cols=DIMSIZE(diam,0)  
////	
////	MAKE/D/N=(rows,cols) cnt_init
////	cnt_init =  cnt
//	//SMOOTH/B= 3 2, cnt
//	SMOOTH/B=2 5,cnt    
////      	
////      	MAKE/O/D/N=(cols) cnt_wave,cnt_fitted_wave        
////	MAKE/O/D/N=(rows,cols) cnt_fitted_matrix
////	FOR(i=0;i<rows;i+=1) 
////      		FOR(j=0;j<cols;j+=1)    
////      			cnt_wave[j] = cnt_init[i][j]	    
////          	ENDFOR       
////          	   
////		DUPLICATE/O/D cnt_wave cnt_wave_dum   
////		CurveFit/NTHR=0 LogNormal  cnt_wave_dum /X=diam  /D=cnt_fitted_wave
////			
////		FOR(j=0;j<cols;j+=1)       
////      			cnt_fitted_matrix[i][j] = cnt_fitted_wave[j]	   
////          	ENDFOR           	  
////	ENDFOR     
////	DUPLICATE/O/D cnt_fitted_matrix cnt   
////	
////	    //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
////	// Note that for dNdlogDp and raw counts matrix should be rows=rows_diam AND cols=rows_t_wave
////	//Since we have opposite formatting as required, we need to transpose dNdlogDp and raw count matrix's
////	 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
////      MAKE/O/D/N=( (DIMSIZE(cnt,1)), (DIMSIZE(cnt,0)) ) cnt_dum
////      DUPLICATE/O/D cnt cnt_dum
////      Redimension cnt_dum
////	MatrixTranspose cnt_dum
////	DUPLICATE/O/D cnt_dum cnt 
////         
////       KILLWAVES cnt_dum 
////       KILLWAVES  X0
//
//END  

FUNCTION Input_dNdlogDp_matrix(ctrlName) : ButtonControl    
  	STRING ctrlName 
	STRING pathname=""  
	STRING filename=""	                                    
	VARIABLE rows, cols,i,j    
  
 	// ################### Load DATE and TIME into a matirx ######################
	LOADWAVE/A=T/J/M/D/A=WAVE/P=$pathname/K=1/V={"\t,"," $",0,1}/L={0,16,0,1,2} filename
 	
 	// MAKE SMPS_datetime wave
	MAKE/O/D/N=(DIMSIZE(T0,0)) SMPS_datetime,t_wave
	WAVE T0,t_wave
	VARIABLE row_index = 0 		
	FOR(i = 0; i < DIMSIZE(T0,0); i+=1)
			SMPS_datetime[i] = T0[row_index][0]  + T0[row_index][1] 
			t_wave[i] = T0[row_index][1] 
			row_index+=1
	ENDFOR
	KILLWAVES/Z  T0  
	//################### READ DIAMETER (Dp) ######################
 	VARIABLE refNum
       filename=INDEXEDFILE($pathname, 0, ".txt")
       //OPEN file FOR read.
	OPEN/R/Z=2/P=$pathName refNum as fileName
	
	//Store results from OPEN in a safe place.
	variable err = V_flag
	STRING Str
	STRING fullPath = S_fileName
	fullPath = S_fileName
	IF (err == -1)
		Print "DemoOPEN canceled by user."
		RETURN -1
	ENDIF
	IF (err != 0) 
		DOAlert 0, "Error in DemoOPEN"  
		RETURN err
	ENDIF
	i = 1
	// Go to the 16th line - this is the last line of header and contains variables names that we need
	FOR(i=1; i<= 16; i+= 1)
		FReadLine refNum, str // Read first line into STRING VARIABLE
	ENDFOR
		VARIABLE start // The starting index of the Bins STRING 
		VARIABLE last // The ENDing index of the Bins STRING 
		VARIABLE result
      start = STRSEARCH(str, "Midpoint", 0)		// Obtaining the starting index of the Bins STRING 
	// Note that start denotes the beginning of the STRING "Midpoint", hence start+9 is the first next character after it
	last = STRSEARCH(str, "Scan Up Time(s)", start+9)  // Obtaining the ENDing index of the Bins STRING 
	
	STRING NumStr_min = str [start+9, last-2]
     
      	result = STRSEARCH(NumStr_min, "\t", Inf, 1) // backwords
	STRING Tail_min = NumStr_min [result+1, STRLEN(NumStr_min)-1]
	      
       //OPEN file FOR read.
	OPEN/R/Z=2/P=$pathName refNum as fileName
	
	// Store results from OPEN in a safe place.
	err = V_flag
	fullPath = S_fileName
	IF (err == -1)
		Print "DemoOPEN canceled by user."
		RETURN -1
	ENDIF
	IF (err != 0)
		DOAlert 0, "Error in DemoOPEN"
		RETURN err
	ENDIF
	i = 1
	// Go to the 16th line
	FOR(i=1; i<= 16; i+= 1)
		FReadLine refNum, str // Read first line into STRING VARIABLE
	ENDFOR
      	start = STRSEARCH(str, "Midpoint", 0)	 	// Obtaining the starting index of the Bins STRING 
	// Note that start denotes the beginning of the STRING "Midpoint", hence start+9 is the first next character after it
	last = STRSEARCH(str, "Scan Up Time(s)", start+9)  // Obtaining the ENDing index of the Bins STRING 
	
	//Storing the subSTRING into NumStr
	STRING NumStr_max = str [start+9, last-2]
       
	//Get the last number in the Num_Str_max
       result = STRSEARCH(NumStr_max, Tail_min, 0) // backwords
	STRING Tail = NumStr_max [result, STRLEN(NumStr_max)-1]
	STRING Head = NumStr_min[0, STRSEARCH(NumStr_min,Tail_min,0)-1]
       STRING Midpoints = Head + Tail
       
       //BREAK the Midpoint STRING into Tockens
       result = 0
       VARIABLE N = STRLEN(Midpoints)
       VARIABLE Count = 0  
       VARIABLE index = 0
	// Determine the length of the entier sting
	DO
		IF((result == -1) || (index >= N)) // The exit condition, when no more cells are found or the index reaches N (END of the STRING)
			Count += 1 
			Print "Count :",  Count 
			BREAK
		ENDIF
		result = STRSEARCH(Midpoints, "\t", index)	// Get the index of "\t"
		IF(result != -1)		// IF found
			Count += 1	// increment Count of cells by 1
			index = result + 1 // set the starting index of next search iteration to the last "\t" index found plus 1 
		ENDIF
	WHILE(1)
	    
	VARIABLE M_Col = Count
	MAKE/O/D/N=(M_Col) Dp
	result = 0
       Count = 0
       index = 0
       i = 0
	// Determine the length of the entier sting
	DO
		IF((result == -1) || (index >= N)) // The exit condition, when no more cells are found or the index reaches N (END of the STRING)
			Count += 1 
			Dp[Count-1] = STR2NUM(Midpoints[index, N-1])
			BREAK
		ENDIF
		result = STRSEARCH(Midpoints, "\t", index)	// Get the index of "\t"
		IF(result != -1)		// IF found
			Count += 1	// increment Count of cells by 1
			Dp[Count-1] = STR2NUM(Midpoints[index, result])
			index = result + 1 // set the starting index of next search iteration to the last "\t" index found plus 1 
		ENDIF
	WHILE(1) 

	//################### READ dNdlogDp matrix  - this is column formatted matrix ######################
   	rows=DIMSIZE(SMPS_datetime,0)
	cols=DIMSIZE(Dp,0)
 	LOADWAVE/A=X/O/J/M/D/A=WAVE/P=$pathName/K=0/L={0,16,0,3,cols} fileName

	//WAVE DC_factor
	//MAKE/O/D/N=(cols) DC_factor_div
	//DC_factor_div = NAN
	//DC_factor_div = DC_factor
      MAKE/O/D/N=(rows, cols) dNdlogDp  
	WAVE X0
      FOR(i=0;i<rows;i+=1) 
      		FOR(j=0;j<cols;j+=1)   
          		dNdlogDp[i][j] = X0[i][j+1]      
          		//IF(DC_factor_div[j])
 	         	//	dNdlogDp[i][j] = dNdlogDp[i][j] / DC_factor_div[j]  
          		//ENDIF
        	ENDFOR 
      ENDFOR         
	
	//######### Reading geometric mean diameter ################
	WAVE Diam
	DUPLICATE/D Dp Diam 
	MAKE/D/N=(NUMPNTS(Diam)) radius 
	radius = Diam/2.
	  
      VARIABLE start_col= 4+count+19  
     LoadWave/A/J/M/D/A=WAVE/P=$pathname/K=1/L={0,16,0,start_col,5} filename
     Duplicate/O/D Wave0, Geo_diam_data      
                        
     MAKE/O/D/N=(DIMSIZE(Geo_diam_data,0) ) Geo_diam, Geo_diam_std
               FOR (i=0;i<DIMSIZE(Geo_diam_data,0) ;i+=1)
	      		Geo_diam[i] = Geo_diam_data[i][2]
	      		Geo_diam_std[i] = Geo_diam_data[i][4]    
		ENDFOR     
      KILLWAVES  X0          
      KILLWAVES wave0        
         
//     DUPLICATE/D Geo_diam Geo_diam_wave  
//     DUPLICATE/D Geo_diam_std Geo_diam_std_wave        
//     SMOOTH/B=3 3, Geo_diam_wave,Geo_diam_std_wave
//
//     MAKE/D/N=(DIMSIZE(Geo_diam,0)) lower_limit,upper_limit,lower_limit_fitted,upper_limit_fitted
//         
//     lower_limit = Geo_diam_wave - (2.0*Geo_diam_std_wave)
//     upper_limit = Geo_diam_wave + (2.0*Geo_diam_std_wave)   
//           
//     CurveFit/M=2/W=0 poly 20, lower_limit/D=lower_limit_fitted    
//     CurveFit/M=2/W=0 poly 20, upper_limit/D=upper_limit_fitted           
//      
//     MAKE/O/D/N=(rows, cols) dNdlogDp_dummy      
//     FOR(i=0;i<rows;i+=1) 
//      		FOR(j=0;j<cols;j+=1)         
//        		
//          		IF(Diam[j]  >= lower_limit_fitted[i] && Diam[j]  <= upper_limit_fitted[i])  
//               		dNdlogDp_dummy[i][j] = dNdlogDp[i][j]
//               	ELSE
//               		dNdlogDp_dummy[i][j] =1E-6
//          		ENDIF   
//          		          
//        	ENDFOR       
//      ENDFOR           
//	DUPLICATE/O/D dNdlogDp_dummy dNdlogDp                 
	   
      	//############# Fitting lognormal fucntion to the data ####################
	SMOOTH/B=2 5,dNdlogDp  
	
	//fit_lognormal()       
	//DUPLICATE/O/D dNdlogDp_fitted_matrix dNdlogDp              
	
	
	//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	// Note that for dNdlogDp and raw counts matrix should be rows=rows_diam AND cols=rows_t_wave
	//Since we have opposite formatting as required, we need to transpose dNdlogDp and raw count matrix's
	 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
      MAKE/O/D/N=( (DIMSIZE(dNdlogDp,1)), (DIMSIZE(dNdlogDp,0)) ) dNdlogDp_dum
      DUPLICATE/O/D dNdlogDp dNdlogDp_dum,dNdlogDp_plot
      Redimension dNdlogDp_dum
	MatrixTranspose dNdlogDp_dum
	DUPLICATE/O/D dNdlogDp_dum dNdlogDp   
	KILLWAVES dNdlogDp_dum   
	
 END  
 

//################### READ Raw Counts matrix  - this is also column formatted matrix ######################
FUNCTION Input_cnt_matrix(ctrlName) : ButtonControl    
  	STRING ctrlName 
	STRING pathname=""  
	STRING filename=""	                                 
	VARIABLE rows, cols,i,j    
	WAVE t_wave,Diam
	WAVE lower_limit_fitted,upper_limit_fitted
	
	rows=DIMSIZE(t_wave,0)
	cols=DIMSIZE(Diam,0)  
    	LOADWAVE/A=X/O/J/M/D/A=WAVE/P=$pathName/K=0/L={0,16,0,3,cols} fileName

      MAKE/O/D/N=(rows, cols) cnt
	WAVE X0
      FOR(i=0;i<rows;i+=1) 
      		FOR(j=0;j<cols;j+=1)  
          		cnt[i][j]=X0[i][j+1]  
        	ENDFOR
      ENDFOR
      	 
    	SMOOTH/B=2 5, cnt        
      	
//      	MAKE/O/D/N=(cols) cnt_wave,cnt_fitted_wave        
//	MAKE/O/D/N=(rows,cols) cnt_fitted_matrix
//	FOR(i=0;i<rows;i+=1) 
//      		FOR(j=0;j<cols;j+=1)  
//      			cnt_wave[j] = cnt_init[i][j]	    
//          	ENDFOR     
//          	   
//		DUPLICATE/O/D cnt_wave cnt_wave_dum   
//		CurveFit/NTHR=0 LogNormal  cnt_wave_dum /X=Dp  /D=cnt_fitted_wave
//		//CurveFit/M=2/W=0 poly 20, cnt_wave_dum /X=Dp  /D=cnt_fitted_wave
//		
//		FOR(j=0;j<cols;j+=1)       
//      			cnt_fitted_matrix[i][j] = cnt_fitted_wave[j]	   
//          	ENDFOR           	  
//	ENDFOR     
//	DUPLICATE/O/D cnt_fitted_matrix cnt   


	
//	MAKE/O/D/N=(rows, cols) cnt_dummy      
//      FOR(i=0;i<rows;i+=1)    
//      		FOR(j=0;j<cols;j+=1)         
//        		
//          		IF(Diam[j]  >= lower_limit_fitted[i] && Diam[j]  <= upper_limit_fitted[i])
//               		cnt_dummy[i][j] = cnt_fitted_matrix[i][j]   
//               	ELSE
//               		cnt_dummy[i][j] =1E-6
//          		ENDIF      
//          		        
//        	ENDFOR       
//      ENDFOR           
//	DUPLICATE/O/D cnt_dummy cnt   

		        
       //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
	// Note that for dNdlogDp and raw counts matrix should be rows=rows_diam AND cols=rows_t_wave
	//Since we have opposite formatting as required, we need to transpose dNdlogDp and raw count matrix's
	 //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
      MAKE/O/D/N=( (DIMSIZE(cnt,1)), (DIMSIZE(cnt,0)) ) cnt_dum
      DUPLICATE/O/D cnt cnt_dum
      Redimension cnt_dum
	MatrixTranspose cnt_dum
	DUPLICATE/O/D cnt_dum cnt 
         
       KILLWAVES cnt_dum 
       KILLWAVES  X0
END       
     


//FUNCTION input_time(ctrlName) : ButtonControl    
//	STRING ctrlName 
//	STRING pathname="" 
//	STRING filename=""	
//	LOADWAVE/A/J/D/A=t_wave/P=$pathname/K=0 filename
//	DUPLICATE/O t_wave0 t_wave  
//	KILLWAVES/Z t_wave0
//END    
//     
//FUNCTION input_diam(ctrlName) : ButtonControl                
//	STRING ctrlName 
//	STRING pathname=""
//	STRING filename=""	
//	LoadWave/A/J/D/A=diam/P=$pathname/K=0 filename 
//	DUPLICATE/O diam0 diam
//	MAKE/D/O/N=(NUMPNTS(diam)) radius
//	radius = diam/2. 
//	KILLWAVES/Z diam0
//END  
// 
// // this input count wave for definition dndlogdp/count
////FUNCTION input_counts(ctrlName) : ButtonControl   
////	STRING ctrlName 
////	STRING pathname=""
////	STRING filename=""	
////	LOADWAVE/J/D/N=cnt/O/K=0/P=$pathname filename 
////	DUPLICATE/O cnt0 cnt
////	KILLWAVES/Z cnt0
////END
//
//// this input count matrix , same as dndlogdp format
//FUNCTION input_counts(ctrlName) : ButtonControl   
//	STRING ctrlName 
//	STRING pathname=""  
//	STRING filename=""	
//	LoadWave/J/M/U={0,0,0,0}/D/O/N=cnt/K=0/P=$pathname filename
//	DUPLICATE/O cnt0 cnt
//	KILLWAVES/Z cnt0
//END  
// 
//
//FUNCTION input_n(ctrlName) : ButtonControl     
//	STRING ctrlName 
//	STRING pathname=""
//	STRING filename=""	
//	LoadWave/J/M/U={0,0,0,0}/D/O/N=dNdlogDp/K=0/P=$pathname filename
//	DUPLICATE/O dNdlogDp0 dNdlogDp
//	KILLWAVES/Z dNdlogDp0  
//END 
//  
 

FUNCTION Makebanana_func(ctrlName) : ButtonControl      
 	STRING ctrlName
 	
 	WAVE t_wave, dNdlogDp_plot, diam,Geo_diam,Geo_diam_std
	Makebanana(t_wave, dNdlogDp_plot, diam,Geo_diam,Geo_diam_std )
	
END
  
//FUNCTION fit_lognormal()   
//	WAVE dNdlogDp,diam,t_wave
// 	VARIABLE rows=DIMSIZE(t_wave,0)   
//	VARIABLE cols=DIMSIZE(diam,0)
//	VARIABLE i,j
//		
//	MAKE/O/D/N=(cols) dNdlogDp_wave,dNdlogDp_fitted_wave
//	MAKE/O/D/N=(rows,cols) dNdlogDp_fitted_matrix
//	FOR(i=0;i<rows;i+=1) 
//      		FOR(j=0;j<cols;j+=1) 
//      			dNdlogDp_wave[j] = dNdlogDp[i][j]	  
//          	ENDFOR 
//		DUPLICATE/O/D dNdlogDp_wave dNdlogDp_wave_dum   
//		
//		CurveFit/NTHR=0 LogNormal  dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
//		//CurveFit/M=2/W=0 poly 20,  dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
//		
//		//CurveFit/NTHR=0 poly 7, dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
//		// CurveFit/NTHR=0 exp  dNdlogDp_wave_dum /X=Dp  /D=dNdlogDp_fitted_wave
////		 display dNdlogDp_wave_dum vs Dp
////		 AppendToGraph dNdlogDp_fitted_wave vs Dp
////		 ModIFyGraph log(bottom)=1 
////		 SetAxis bottom 3.2,120	  
////		 ModifyGraph mode=4,rgb=(0,0,0) 	  
////		PauseUpdate    
//		FOR(j=0;j<cols;j+=1)    
//      			dNdlogDp_fitted_matrix[i][j] = dNdlogDp_fitted_wave[j]	   
//          	ENDFOR 
//          	  
//	ENDFOR       
//	        
//END          
//     

  
FUNCTION Makebanana(t_wave, dNdlogDp_plot, diam,Geo_diam,Geo_diam_std )
    	WAVE t_wave, dNdlogDp_plot, diam ,Geo_diam,Geo_diam_std
    	DUPLICATE/O t_wave time_plotw
    	DUPLICATE/O diam diam_plotw
	DUPLICATE/O dNdlogDp_plot n_plotw 
    	n_plotw = n_plotw/1000.
    	
    	VARIABLE N_Rows = DIMSIZE(time_plotw,0)	
	VARIABLE N_Cols = DIMSIZE(n_plotw,1)	
		 
   	Display;AppENDMatrixContour n_plotw vs {time_plotw, diam_plotw}
	Removecontour n_plotw   
	REDIMENSION/D/N=(N_Cols+1) diam_plotw 
	REDIMENSION/D/N=(N_Rows+1) time_plotw	
	diam_plotw[N_Cols] = diam_plotw[N_Cols-1] + 5
	time_plotw[N_Rows] = time_plotw[N_Rows-1] + 10
	AppENDImage n_plotw vs {time_plotw, diam_plotw}
	Setscale d, 0,1,time_plotw
	ModIFyGraph margin(right) = 90    
	ModIFyGraph nticks=5
	ModIFyGraph log(left)=1 
	ModifyImage n_plotw ctab= {0,15,Rainbow,1}
	ColorScale/C/N=text0/F=0/A=MC/X=68/Y=0 HeightPct=116
	ColorScale/C/N=text0"\F'Times New Roman'\Z12dN/dlogDp (10\S3\M\Z12 cm\S-3\M\Z12)"
	//label bottom "Time (EST)"
	Label left "\F'Times New Roman'\Z12Dp (nm)"
	ModIFyGraph nticks(left)=5  	 
	//AppendToGraph Geo_diam vs t_wave
	AppendToGraph/T Geo_diam vs t_wave
	ModifyGraph mode=4,rgb=(0,0,0)
	//ModifyGraph mode=4,rgb=(0,0,0)
END      
   

FUNCTION Select_NPFtime(ctrlName) : ButtonControl    
  	STRING ctrlName 
  	
 	WAVE t_wave, dNdlogDp, diam, cnt
	NVAR NPFSTime, NPFEtime
      VARIABLE NPF_num_time = NPFEtime - NPFStime+1   
	PRINT NPF_num_time  
	  
	MAKE/O/D/N=(NPF_num_time) NPF_time
	MAKE/O/D/N=(DIMSIZE(diam,0),NPF_num_time) NPF_dNdlogDp , NPF_cnt
	VARIABLE i,j
	VARIABLE k=0
	FOR(i=NPFStime;i<=NPFEtime;i+=1)
  		NPF_time[k]=t_wave[i]
		k+=1          
	ENDFOR    
   
	k=0  
	FOR(i=0;i<DIMSIZE(diam,0);i+=1)
		FOR(j=NPFStime;j<=NPFEtime;j+=1)
   			NPF_dNdlogDp[i][k] = dNdlogDp[i][j]     
          		NPF_cnt[i][k] = cnt[i][j] 
         	  	k+=1      
		ENDFOR   
	       k=0      
   	ENDFOR    
   	  
   	DUPLICATE/O/D NPF_cnt cnt 
   	DUPLICATE/O/D NPF_time t_wave
   	DUPLICATE/O/D  NPF_dNdlogDp dNdlogDp
   	Save/J/M="\r\n" NPF_time as "NPF_time.txt"   
END   
   

//####################################################################################
//############################ Obtaining Growth Rates #####################################
//#####################################################################################
 
 
FUNCTION Growth_Rate(ctrlName) : ButtonControl      
 	STRING ctrlName
	calc_r_bound()
	starting()
	 turn_size_order_around()
	 make_init_waves()
	cal_losses()
	Call_Fitting_dNc_dt()
	Call_Fit_range()
END
 
   
END

//FUNCTION calc_r_bound(ctrlName) : ButtonControl      
 	//STRING ctrlName
FUNCTION calc_r_bound()
 	WAVE radius
	VARIABLE d_logr_avg
	MAKE/D/O/N=(NUMPNTS(radius)+1) r_bound
	r_bound = 10^( (LOG(radius[p]) + LOG(radius[p+1]) ) /2 )
	DUPLICATE/O/D radius d_logr 
	d_logr = LOG(r_bound[p+1]) - LOG(r_bound[p])
	d_logr_avg = Mean(d_logr, 1, (NUMPNTS(d_logr)-1) )

	MAKE/D/O/N=(NUMPNTS(radius)+1) r_bound
	DUPLICATE/O/D radius d_logr
	r_bound[0] = 10^(LOG(r_bound[1]) - d_logr_avg )
	r_bound[(NUMPNTS(r_bound)-1)] = 10^(LOG(r_bound[(NUMPNTS(r_bound)-2)]) + d_logr_avg ) 
	d_logr = LOG(r_bound[p+1]) - LOG(r_bound[p])
	d_logr_avg = Mean(d_logr, 1, (NUMPNTS(d_logr)-1) )
END
	 
	
//FUNCTION starting (ctrlName) : ButtonControl
	//STRING ctrlName
FUNCTION starting()
	VARIABLE scan 
	start(0)    
END


FUNCTION start (scan)     
	VARIABLE scan
	VARIABLE/G inc 
	inc=1
	// Function 'start' creates the input file for the calculations, taken from the n- and t-matrix.
	// Dependent on the lay-out of the input files, changes may be necessary.
	// Units can be converted in this macro, or before entering, as preferred. 
	// If size bins are going from large to small sizes, choose "bin=p"
	// If size bins are going from small to large sizes, choose "bin=L-1-p"
	// This will create a wave with bin numbers that increase with increasing size (starting at bin=0)
	// Bin size needs to be expressed as radius in nanometers. 
	// Conversion factor can be applied where waves "radius" and "r_bound" are being read in. 
	// Global variables (which are accessed regularly throughout the program) are defined 
	// at start of this macro only, their values should only be changed here.
	VARIABLE/G Temp				// Temperature
	VARIABLE/G density				// density of particle
	VARIABLE/G pressure				// pressure 
	VARIABLE/G Boltz				// Boltzmann's constant
	VARIABLE/G Kwconst				// constant for in wall loss rate expression
	VARIABLE/G speed				// mean molecular speed of air
	VARIABLE/G speed_gas			// mean molecular speed of condensing species
	VARIABLE/G diff_gas				// gas-phase diffusion coefficient of condensing species
	VARIABLE/G S_V_ratio			// [m^-1] surface-to-volume ratio in chamber OR inverse of mixed layer depth in atmosphere
	VARIABLE/G f_s					// ratio of projected horizontal surface area to total surface area in chamber
	VARIABLE/G G					// acceleration of gravity
	VARIABLE/G r_nucl				// assumed radius of critical cluster
	VARIABLE/G MW_H2SO4			// molecular weight of H2SO4
	VARIABLE/G MW_air				// molecular weight of air
	VARIABLE/G MW_organic			// molecular weight of organic condensable species (i.e. pinox)
	VARIABLE/G gas_const			// molar gas constant
	VARIABLE/G MW_organic			// molecular weight of organic condensable species (i.e. pinox)
	VARIABLE/G gas_const			// molar gas constant
	VARIABLE/G alpha				// accommodation coefficient
	
	VARIABLE/G L 	 				// number of size bins in each input wave (i.e. number of rows in input matrix)
	VARIABLE/G K	 				// number of measurement times (i.e. number of columns in input matrix)
	VARIABLE/G r_det				// minimum detectable radius
	WAVE radius, r_bound, dNdlogDp, t_wave	// these waves should be imported.
	WAVE wsu_t
	//wave sd_n, cnt			 		// one of these two waves should be imported; other does not exist.
	//wave n_back_all, r_back_all 		// these waves only exist after analysis has been performed.
	//wave dN_dlnr1_test_recent_nucl, dN_dlnr2_test_recent_nucl
	//wave n1_test_1nm_wide, n1_test_1nm_noise, n1_test_1nm_cutoff_wide, n1_test_1nm_cutoff_noise
	//wave n2_test_1nm_wide, n2_test_1nm_noise, n2_test_1nm_cutoff_wide, n2_test_1nm_cutoff_noise
	
	
	//DUPLICATE/O/D synch_t_wsu wsu_t
	//Interpolate2/T=2/N=(K)/Y=wsu_t synch_t_wsu

	Temp = 295.25  		// wsu_t[scan] + 273  //295.0			// Kelvin		
	density = 1.76			// g/cc (stimated at 1.0 for organic aerosol; 1.2 for inorganic aerosol)	  	
	pressure = 1000.00		// mbar
	Boltz = 1.38e-16		// erg/K = 1.38e-23 J/K (erg = g cm^2 s^-2)	
	Kwconst = 3.6e5 		// 1e4 // atmospheric guess                 
	MW_air = 28.97	 	// g/mol  	
	MW_H2SO4 = 98		// 132.76 for (NH4)2SO4      // g/mol
	MW_organic = 100	// g/mol (estimate)
	S_V_ratio =  0.000017  // i.e. 1/1000 m	= 0.001 m-1 =1E-3 m-1 = 1E-5 cm-2 cm-1 for mixed layer depth of 1000m
	f_s = 1.0				// unitless. Atmosphere=1  Calspan=0.167  PSI=0.167 I think (1/6)
	G = 980.0			// [cm/s^2]
	gas_const = 8.314		// J mol^-1 K^-1
	alpha = 1.0			// unitless
	 
	speed = SQRT(8*8.314E7*Temp/(pi*MW_air))				// cm/s	for AIR only     `
	speed_gas = SQRT(8*8.314E7*Temp/(pi*MW_H2SO4))		// cm/s	for H2SO4
 
	// note that RADIUS is used, not diameter!! Convert if necessary.
	// Likewise, n = dN/dln(R) , nn = dN/dR
	// Conversions: dN/dln(R) = dN/dR * R = dN/dlog(R) / 2.303 = dN/dln(D) = dN/dD * D 
	//DUPLICATE/O/D n nn    
	//nn = n/radius[p]  
	r_nucl = 0.5 				// assumed radius of critical cluster
	r_det = radius[L-1]			// minimum detectable radius
	L=NUMPNTS(radius)	 	// number of size bins (i.e. number of rows in input matrix)
	K = DIMSIZE(dNdlogDp,1)			// number of measurement times (i.e. number of columns in input matrix)
	MAKE/N=(L) /D/O bin		//creates a wave of size equal to number of size bins 
	bin = p					// sets values of bin equal to the point number [p], from small to large values; 
	DUPLICATE/O/D bin n_t1, n_t2, Nc_t1, Nc_t2, nn_t1, nn_t2, n_avg, Nc_avg, nn_avg
	DUPLICATE/O/D bin ln_r, d_r, d_lnr, ln_r, dn_dt, dNc_dt, dn_dr_t1, dn_dr_t2, dn_dr_avg
	DUPLICATE/O/D dNdlogDp V_dist
	DUPLICATE/O/D bin V_dist_t1, V_dist_t2, dV_dlnr1, dV_dlnr2
	//Make/O/D/N=(K*L) n_back_t1, n_back_t2, r_back_t1, r_back_t2
	MAKE/N=1/D/O dt			// seconds
	//Duplicate/O/D n_t1 n_uncorr_t1, n_uncorr_t2 
	ln_r = ln(radius) 
	d_r = r_bound[p]-r_bound[p+1]
	d_lnr = ln(r_bound[p])-ln(r_bound[p+1])			// d(ln(R)) using the bin boundaries	
	V_dist = dNdlogDp[p][q]*d_lnr[p]*(4/3)*PI*(radius[p])^3*1e-9	// Volume in each size bin (um^3)
	V_dist_t1 = V_dist[p][scan]
	V_dist_t2 = V_dist[p][scan+inc]
	dV_dlnr1 = V_dist_t1 / d_lnr
	dV_dlnr2 = V_dist_t2 / d_lnr
	n_t1 = dNdlogDp[p][scan]		// to read in from dN/dlnr matrix (produced by awk)
	n_t2 = dNdlogDp[p][scan+inc]	// dN/dlnr (units of r don't matter)

	replace_NaN_by_num(n_t1, 0) 
	replace_NaN_by_num(n_t2, 0)
	 
	//Duplicate/O n_t1,n_t1_smth
	SMOOTH/B=2 5, n_t1
	//Duplicate/O n_t2,n_t2_smth  
	SMOOTH/B=2 5, n_t2
	
	
	
	nn_t1 = n_t1/radius	// nn = dN/dR
	nn_t2 = n_t2/radius
	//n_uncorr_t1 = n_uncorrected[p][scan]
	//n_uncorr_t2 = n_uncorrected[p][scan+1]
	//n_back_t1 = n_back_all[p][scan] 
	//n_back_t2 = n_back_all[p][scan+1]
	//r_back_t1 = r_back_all[p][scan]
	//r_back_t2 = r_back_all[p][scan+1]
	// procedure for rectangular integration of XY data
	MAKE/D/O/N=(L) n_d_lnr_1, n_d_lnr_2, R_ip, delta_r
	//n_d_lnr_1 = n_t1_smth*d_lnr			// product of integration
	n_d_lnr_1 = n_t1*d_lnr					// product of integration
	Nc_t1 = SUM(n_d_lnr_1, 0, p)			// summation cq. rectangular integration
	//n_d_lnr_2 = n_t2_smth*d_lnr			// The imported Nc is a straight summation, so this method...
	n_d_lnr_2 = n_t2*d_lnr					// The imported Nc is a straight summation, so this method...
	Nc_t2 = SUM(n_d_lnr_2, 0, p)			// ...of integration is more suitable then trapezoidal.
	dt = t_wave[scan+inc] - t_wave[scan]		// time interval (seconds) between two consecutive measurements
	make/D/O/N =(K-1) t_avg_wave			// end times of size dist measurement
	t_avg_wave = (t_wave[p]+t_wave[p+inc]) / 2
	duplicate/O t_wave t_midnucl
	duplicate/O t_avg_wave t_midgrowth
	t_midnucl = t_wave  +115
	t_midgrowth = (t_midnucl[p]+t_midnucl[p+inc]) / 2 
	 
	//K=NUMPNTS(t_wave)	  
	//variable i
	//MAKE/D/O/N =(K-1) t_avg_wave			// end times of size dist measurement
	//FOR(i=0;i<K;i+=1)
	//	t_avg_wave[i] = (t_wave[i]+t_wave[i+inc]) / 2
	//ENDFOR    
	//DUPLICATE/O t_wave t_midnucl
	//t_midnucl = t_wave + 40
	//MAKE/D/O/N =(K-1) t_midgrowth 
	//FOR(i=0;i<K;i+=1)
	//	t_midgrowth[i] = (t_midnucl[i]+t_midnucl[i+inc]) / 2 
	//ENDFOR
	   
	dn_dt = (n_t2 - n_t1) / dt	   				// change in (nat.) logarithmic particle number distribution function in time
	//n_avg = (n_t1_smth + n_t2_smth) / 2
	n_avg = (n_t1 + n_t2) / 2
	dNc_dt = (Nc_t2 - Nc_t1) / dt		 		// change of cumulative number in time: this is the quantity being fit
	Nc_avg = (Nc_t1 + Nc_t2) / 2
	nn_avg = (nn_t1 + nn_t2) / 2
	R_ip = INTERP(Nc_t2, Nc_t1, radius)		// Using interpolation to obtain estimate of g = delta_r / dt
	delta_r = radius-R_ip	 	 			// the change in radius over time for each cohort of particles
	differentiateXY (radius, n_t1, "dn_dr_t1") 	// central difference
	differentiateXY (radius, n_t2, "dn_dr_t2")                               
	//dn_dr_t1 = (n_t1[p+1]-n_t1[p]) / (radius[p+1]-radius[p]) // backwards difference
	//dn_dr_t2 = (n_t2[p+1]-n_t2[p]) / (radius[p+1]-radius[p])
	//dn_dr_t1[0] = 0		 // because backwards difference gives NaN which could screw things up.
	//dn_dr_t2[0] = 0
	dn_dr_avg = (dn_dr_t1 + dn_dr_t2) / 2
	 
	//stdev_sd_n(scan)   
	stdev_cnt(scan)  
	//stdev_prop(scan)      
	
	//KILLWAVES n_d_lnr_1, n_d_lnr_2, dn_dr_t1, dn_dr_t2//, V_dist_t1, V_dist_t2 
	
	replace_NaN_by_num(Nc_t2, 0) 
	replace_NaN_by_num(Nc_t1, 0)
	
	MAKE/O/D/N=(k) delta_N_by_N
	MAKE/O/D/N=(k-1) delta_n_dum
	
	delta_N_by_N = ( ( SUM(Nc_t2, 0, p) - SUM(Nc_t1, 0, p)  ) / SUM(Nc_t1, 0, p) )
	//delta_n_dum = delta_N_by_N
	//print scan, delta_N_by_N[scan+inc]
	delta_n_dum[scan] = delta_N_by_N[scan+inc]
	//print delta_n_dum[scan]
	//print "dum", delta_n_dum	 
END        
  
 

FUNCTION replace_nan_by_num(name_of_wave, num)	// replaces -999 or less (as output by labview averaging program) by NaN.
	VARIABLE num								// number to replace NaN by.
	WAVE name_of_wave
	VARIABLE i
	i = 0
	DO											//if (name_of_wave[i] > -10) // i.e. if it is a real data value		
		IF (NUMTYPE(name_of_wave[i]) == 0) 		// i.e. if it is a real data value
			name_of_wave[i] = name_of_wave[i]
		ELSE 									// i.e. if it is either INF or NaN
			name_of_wave[i] = num
		ENDIF
		i = i+1
	WHILE(i < NUMPNTS(name_of_wave))	
END


FUNCTION stdev_cnt(scan) 			// Calculates sd(dNc/dt) from imported matrix of # counts
	VARIABLE scan
	WAVE bin, cnt, Nc_t1, Nc_t2, dNc_dt, dt, n_t1, d_lnr, Nc_avg, dNdlogDp, n_t2  //, prop_11_16_run13
	NVAR L 	
	VARIABLE f_sd, i
	DUPLICATE/O/D bin cnt_t1, cnt_t2, prop_t1, prop_t2, var_dn, var_dNc, sd_dNc_dt, sd_n
	DUPLICATE/D/O bin var_stat_t1, var_stat_t2, var_flow_t1, var_flow_t2, var_one, sd_dn_dt, var_dN_ind
	f_sd = 0.01 					// empirical number
	cnt_t1 = cnt[p][scan] 			// number of counts in a bin: taken from matrix of counts as input.
	cnt_t2 = cnt[p][scan+1]
	i = 0
	IF(NUMTYPE(cnt_t1[p]) == 0 )
		DO
			IF ( cnt_t1[i]>0 ) 
				prop_t1[i] = n_t1[i]*d_lnr[i] / cnt_t1[i]  
			ELSE				// proportionality constant between concentration and counts
				prop_t1[i] = 0 		// because otherwise st.dev. becomes NaN.
			ENDIF
			
			IF(cnt_t2[i]>0 )
				prop_t2[i] = n_t2[i]*d_lnr[i] / cnt_t2[i]
			ELSE
				prop_t2[i] = 0
			ENDIF 
			i = i+1  
		WHILE(i<L)
		var_stat_t1 = prop_t1^2*cnt_t1
		var_stat_t2 = prop_t2^2*cnt_t2
	ELSE  //if ( numtype(cnt_t1[p]) == 2 )   // print "a wave 'prop_n' must be created, where prop_n = particle number density / counts"
		//WAVE prop_n 					// print "then activate the 5 lines in Function stdev() dealing with prop_n"
		//prop_t1 = prop_n				// print "and execute same function again"
		//prop_t2 = prop_n
		//var_stat_t1 = prop_n * n_t1*d_lnr
		//var_stat_t2 = prop_n * n_t2*d_lnr
	ENDIF	
	var_one = (14*(prop_t1+prop_t2)/2)^2
	var_flow_t1 = (f_sd*n_t1*d_lnr)^2 
	var_flow_t2 = (f_sd*n_t2*d_lnr)^2
	var_dN = var_stat_t1 + var_stat_t2 + var_flow_t1 + var_flow_t2 + 2*var_one
	var_dN_ind = var_stat_t1 + var_flow_t1 + var_one
	var_dNc = SUM(var_dN, 0, p)
	sd_dNc_dt = SQRT(var_dNc) / dt
	sd_dn_dt = SQRT(var_dN) / (dt*d_lnr)
	KILLWAVES var_stat_t1, var_stat_t2, var_flow_t1, var_flow_t2, var_one, prop_t1, prop_t2
	KILLWAVES cnt_t1, cnt_t2, var_dn, var_dNc
	// the abs. st. dev. is used as weights in fitting procedure and as error bars around dNc/dt in graph.
END

FUNCTION stdev_prop(scan) 	
// Calculates sd(dNc/dt) from imported proportionality wave
VARIABLE scan
WAVE bin, Nc_t1, Nc_t2, dNc_dt, dt, n_t1, d_lnr, Nc_avg, dNdlogDp, n_t2, prop_n,cnt //11_16_run13
NVAR L //number
VARIABLE f_sd, i
DUPLICATE/O/D cnt prop_n
Duplicate/O/D bin, var_dn, var_dNc, sd_dNc_dt, sd_n
duplicate/D/O bin var_stat_t1, var_stat_t2, var_flow_t1, var_flow_t2, var_one, sd_dn_dt, var_dN_ind
f_sd = 0.01 // empirical number
i = 0 
//	print "wave of counts is NOT available" 
//	print "a wave 'prop_n' must be created, where prop_n = particle number density / counts"
//	print "then activate the 5 lines in Function stdev() dealing with prop_n"
//	print "and execute same function again"
	var_stat_t1 = prop_n * n_t1*d_lnr
	var_stat_t2 = prop_n * n_t2*d_lnr
	var_one = 1*(prop_n)^2
	var_flow_t1 = (f_sd*n_t1*d_lnr)^2
	var_flow_t2 = (f_sd*n_t2*d_lnr)^2
	var_dN = var_stat_t1 + var_stat_t2 + var_flow_t1 + var_flow_t2 + 2*var_one
	var_dN_ind = var_stat_t1 + var_flow_t1 + var_one
	var_dNc = sum(var_dN, 0, p)
	sd_dNc_dt = sqrt(var_dNc) / dt
	sd_dn_dt = sqrt(var_dN) / (dt*d_lnr)
	KILLWAVES var_stat_t1, var_stat_t2, var_flow_t1, var_flow_t2, var_one
	KILLWAVES var_dn, var_dNc
// the abs. st. dev. is used as weights in fitting procedure 
// and as error bars around dNc/dt in graph.
END

FUNCTION turn_size_order_around()  //changes order of size around for matrix n and waves radius, r_bound, cnt and/or sd_n (i.e. all input waves)
//FUNCTION turn_size_order_around(ctrlName) : ButtonControl     
 	//STRING ctrlName 	
	VARIABLE/G L
	WAVE radius,r_bound,dNdlogDp,sd_dNc_dt,sd_dn_dt,cnt
	L = NUMPNTS(radius)
	DUPLICATE/O radius r_opp
	DUPLICATE/O r_bound rb_opp
	DUPLICATE/O dNdlogDp n_opp
	DUPLICATE/O sd_dNc_dt sd_dNc_dt_opp
	DUPLICATE/O sd_dn_dt sd_dn_dt_opp
	DUPLICATE/O cnt cnt_opp
	r_opp = radius[L-1-p]
	rb_opp = r_bound[L-p]
	n_opp = dNdlogDp[L-1-p][q] 
	sd_dNc_dt_opp = sd_dNc_dt[L-1-p] 
	sd_dn_dt_opp = sd_dn_dt[L-1-p]
	cnt_opp = cnt[L-1-p][q]

	DUPLICATE/O r_opp radius
	DUPLICATE/O rb_opp r_bound
	DUPLICATE/O n_opp dNdlogDp
	DUPLICATE/O sd_dNc_dt_opp sd_dNc_dt
	DUPLICATE/O sd_dn_dt_opp sd_dn_dt
	DUPLICATE/O cnt_opp cnt	
	KILLWAVES/Z r_opp, rb_opp,n_opp, sd_dNc_dt_opp,sd_dn_dt_opp ,cnt_opp
END 


//FUNCTION make_init_waves(ctrlName) : ButtonControl     
 	//STRING ctrlName 	
FUNCTION make_init_waves()
	MAKE/D/O/N=(8) w
	MAKE/D/O/T/N=(9) fit_param_txt   
	NVAR K
	fit_param_txt[0] = "g over gamma"  
	fit_param_txt[1] = "alpha" 
	fit_param_txt[2] = "diffusion coeff"
	fit_param_txt[3] = "coag all"
	fit_param_txt[4] = "coag loss" 
	fit_param_txt[5] = "coag prod"  
	fit_param_txt[6] = "wall loss diff"
	fit_param_txt[7] = "wall loss grav"
	fit_param_txt[8] = "chi squared"   
	w[0] = 4.6 	//4.6	// 5 
	w[1] = 1  
	w[2] = 0.1 
	w[3] = 1   
	w[4] = 1  
	w[5] = 1
	w[6] = 0         
	w[7] = 0    
	w[8] = 1
	MAKE/D/O/N=((NUMPNTS(w)+1),(K-1)) w_all          
	DUPLICATE/D/O dNc_dt dNc_dt_calc, rel_res, res_dNc_dt
END  
 
  
//FUNCTION  cal_losses(ctrlName) : ButtonControl 
	//STRING ctrlName 
FUNCTION  cal_losses()
	losses()      
END   

  
FUNCTION losses()     
	WAVE radius, bin, n_avg        
	// returns auxiliary quantities necessary to calculate processes affecting dn/dt
	DUPLICATE/O bin B, D_B, K_wall_diff, K_wall_grav, wall_diff_rate, wall_grav_rate
	DUPLICATE/O bin coag_loss_rate, coag_prod_rate //, int_Kc_n//, wall_diff_rate_d  
	Duplicate/O/D bin dn_dt_tubing_loss
	B = mobility(radius)					// mobility (s/g)
	D_B = diffusioncoeff(radius)				// Brownian diffusion coefficient (cm^2/s)
	K_wall_diff = wall_loss_diff(radius)		// diffusional wall loss rate constant ( /s)
	K_wall_grav = wall_loss_grav(radius)		// gravitational wall loss rate constant ( /s)
	wall_loss_diff_rate()
	wall_loss_grav_rate()
	coagulation_loss_rate()
	coagulation_prod_rate()
	first_order_rate() 
	dn_dt_tubing_loss = 4*n_avg / ( 3*(radius)^(7/3) ) // (Eq 39.2 Fuchs)
END 
  

FUNCTIOn mobility(R)    //calculates mobility according to Mike and thus probably Fuchs
	VARIABLE R										// nm
	VARIABLE Viscosity, Lambda, Knudsen, Slip, rho_c
	NVAR Temp, pressure
	//VARIABLE/G B_kin_const
	rho_c = 1000*pressure*SQRT(8*28.97/(PI*8.314e7*Temp))					  // g/(s cm^2)
	Viscosity = (1.718 + 0.0049*(Temp-273.15))*1.0e-4							  // poise (g/cm s)
	Lambda = 1.0e4*Viscosity / (0.499*pressure*SQRT(8*28.97/(PI*8.314e7*Temp))) // nm
	Knudsen = Lambda / R													  // dimensionless
	Slip = 1 + 1.246*Knudsen + 0.42*Knudsen*EXP(-0.87/Knudsen)
	// calculate the kinetic limit mobility: B(kin)=B_kin_const/R^2 					  // (R in nm)
	// for use in backwards calc to integrate Kc(a)da. Checked 28-05-02.
	//B_kin_const = (0.5553*lambda*1e-7) / (2*PI*Viscosity*(1e-7)^2)
	RETURN Slip/(6*PI*Viscosity*R*1e-7)									  // mobility s/g
END

FUNCTION diffusioncoeff(R)	   // Diffusion coefficient of a particle with radius R.
	VARIABLE R
	NVAR Temp, Boltz
	RETURN Boltz * Temp * mobility(R)	// Brownian Diffusion coeffcient cm^2/s
	
END


FUNCTION wall_loss_diff(R)		// first order rate constant for diffusional wall losses
	VARIABLE R					
	VARIABLE rho_c
	NVAR Boltz, pressure, Temp, speed, Kwconst	   
	// temperature dependence, as well as dependence on chamber (or PBL) dimensions, based on Hoppel et al. 2001
	VARIABLE/G Kw_kin_const
	rho_c = 1000*pressure*SQRT(8*28.97/(PI*8.314e7*Temp))		// g/(s cm^2)
	Kw_kin_const = Kwconst*SQRT(Boltz*Temp/rho_c)	 			// cm^-1 s^-0.5
	//Kw_kin_const = 3.6E-3
	
	RETURN SQRT(diffusioncoeff(R)) * Kw_kin_const				// particle (s^-1)
	//return sqrt(diffusioncoeff(R)) * Kwconst * sqrt(Boltz*Temp/rho_c)	// particle (s^-1)
	//return sqrt(0.1) * Kw_kin_const *w[6]							// H2SO4 (g) (s^-1)
END
    

FUNCTION wall_loss_grav(R)	// first order rate constant for gravitational wall losses
	VARIABLE R
	NVAR S_V_ratio
	NVAR f_s  
	NVAR G
	RETURN mass_particle(R) * mobility(R) * G * f_s * S_V_ratio 	// s^-1  //1E-2 is to convert S_V_ratio m-1 to cm-1
END


FUNCTION wall_loss_diff_rate()
	NVAR L
	WAVE ln_r, K_wall_diff, n_avg, wall_diff_rate, d_lnr
	MAKE/O/D/N=(NUMPNTS(ln_r)) wall_int_diff, Kw_n_diff
	Kw_n_diff = K_wall_diff * n_avg
	//IntegrateXYPro (ln_r, Kw_n_diff, "wall_int_diff")
	//wall_diff_rate = abs(wall_int_diff)
	wall_int_diff = Kw_n_diff * d_lnr
	wall_diff_rate = SUM( wall_int_diff, 0, p )
END


FUNCTION wall_loss_grav_rate()
	NVAR L
	WAVE ln_r, K_wall_grav, n_avg, wall_grav_rate, d_lnr
	MAKE/O/D/N=(NUMPNTS(ln_r)) wall_int_grav, Kw_n_grav
	Kw_n_grav = K_wall_grav * n_avg
	//IntegrateXYPro (ln_r, Kw_n_grav, "wall_int_grav")
	//wall_grav_rate = abs(wall_int_grav)
	wall_int_grav = Kw_n_grav * d_lnr
	wall_grav_rate = SUM( wall_int_grav, 0, p )
	KILLWAVES wall_int_grav 	//, Kw_n_grav
END


FUNCTION coagulation_loss_rate()
	NVAR L
	// calculates the loss rate by coagulation "coag_loss_int"
	// -2 Integral( n(a) * Integral( Kc*n(a') ) dlna' ) dlna
	//int_Kc_n is the output of function "coagulation_loss"
	// creates a wave of the loss rate: the outer integral, by applying it to all values of bin e.a.
	WAVE n_avg, coag_loss_rate, bin, d_lnr
	MAKE/O/D/N=(NUMPNTS(d_lnr)) coag_loss_int, integrand_c_l, dn_dt_coag_loss
	integrand_c_l = n_avg * coagulation_loss(bin)
	dn_dt_coag_loss = integrand_c_l * 2 				// dn/dt by coag loss
	//IntegrateXYPro (ln_r, integrand_c_l, "coag_loss_int")
	//coag_loss_rate = 2*abs(coag_loss_int)
	coag_loss_int = integrand_c_l * d_lnr
	coag_loss_rate = 2 * SUM( coag_loss_int, 0, p ) 	// dNc/dt by coag loss
	KILLWAVES coag_loss_int, integrand_c_l
END


FUNCTION coagulation_loss(bin)  	
	// calculates Integral( Kc*n(a') ) dlna'  to be used in function "coagulation_loss_rate"
	// int_Kc_n = Integral( Kc*n(a') ) dlna'  ("the inner integral") 
	VARIABLE bin  			// returns the value of the integral for a particular value of bin. 
	WAVE n_avg, radius, d_lnr
	DUPLICATE/O/D radius Kc_n, int_Kc_n_wave, int_Kc_n_sum	// Kc_n is the integrand
	VARIABLE a1, int_Kc_n
	NVAR L 						// Total number of size bins (30 for Pacific DMA)
	a1 = radius[bin]				// loss of particles of radius a1 being calculated.
	// For a constant value of radius a1, as defined by the bin number, and a changing value of a2, as defined by 
	// looping through radius wave, a wave is created of the product of Kc(a1,a2)*n(a2), called Kc_n
	// Since radius and n_avg are waves, Kc_n will be a wave,  while K-coag is called for constant radius a1 as defined by bin number. 
	Kc_n = coag_rate_constant(a1, radius) * n_avg
	// The integration returns a wave of integrals of Kc_n above a certain value of radius.
	//IntegrateXYPro (ln_r, Kc_n, "int_Kc_n_wave").  We need the total integral of Kc_n for all values of radius, 
	// which is the last value in the wave (=integral above the smallest radius).  Its absolute value is returned to correct the sign (integral must be positive).
	//int_Kc_n = abs(int_Kc_n_wave[L-1])
	int_Kc_n_wave = Kc_n * d_lnr		// what is being integrated
	int_Kc_n_sum = SUM( int_Kc_n_wave, 0, p )
	int_Kc_n = ABS(int_Kc_n_sum[L-1])			
	RETURN int_Kc_n
	//KILLWAVES int_Kc_n_wave, int_Kc_n_sum, Kc_n
END


FUNCTION coagulation_prod_rate()
	VARIABLE i, j
	WAVE bin, radius, coag_prod_rate, d_lnr	//, dn_dt_ind
	NVAR L
	DUPLICATE/D/O radius dn_dt_a, dn_dt_coag_prod, coag_prod_int, dn_dt_b, dn_dt_c
	MAKE/O/D/N=(L,L) dn_dt_matrix
	dn_dt_matrix = 0
	i=0
	DO
		DUPLICATE/D/O radius dn_dt_ind
		coagulation_prod(i)
		dn_dt_matrix[][i] = dn_dt_ind[p]
		// matrix of individual contributions to dn/dt
		// column gives dn/dt for constant reacting particle size a1,  row gives dn/dt for constant formed particle size a
		i = i+1
	WHILE(i<L)
	j=0
	DO
		// dn_dt_a = dn_dt_matrix[j][p] *d_lnr // case 1; d_lnr here means integrating over d_lnr1
		dn_dt_a = dn_dt_matrix[j][p] 					// cases 1*, 2 & 3
		dn_dt_b = dn_dt_matrix[j][p]					// TEST only
	  	// writes row of matrix to wave: dn/dt forming particle size a
		dn_dt_coag_prod[j] = SUM(dn_dt_a, 0, L) 		// dn(a)/dt * d_lnr1 for case 1; dn(a)/dt for case 2 & 3
		dn_dt_c[j] = SUM(dn_dt_b, 0, L)				// TEST only
		// summing each wave gives total dn/dt forming particle size a 
		j = j+1
	WHILE(j<L)
	//coag_prod_int = dn_dt_coag_prod 				// case 1
	coag_prod_int = dn_dt_coag_prod * d_lnr 			// case 2 & 3; d_lnr here means integrating over dln(a)
	coag_prod_rate = SUM(coag_prod_int, 0, p) 		// dNc/dt by coag production, for both cases 1 & 2.
	//KILLWAVES dn_dt_a, coag_prod_int, dn_dt_b, dn_dt_c, dn_dt_matrix
END



FUNCTION coagulation_prod(bin)
	// Calculates the production of particles due to coagulation of a particles with a fixed size, a1 (as defined by bin), with a particle of a range of 
	// sizes, a2 (as defined by radius). Both a1 and a2 are equal to bin radii, and resulting particles of size a are divided over two neighbouring bins. 
	// Note that this does NOT calculate the production of particles of size defined by bin. (Exponential) bin sizes need not be linearly spaced.
	VARIABLE bin
	WAVE radius, n_avg, r_bound, d_lnr, dn_dt_ind //, r_back_ind
	NVAR K, L, r_var
	DUPLICATE/O/D radius a2, a, n_avg_a2, Kc_n_n_J, int_Kc_n_n_J_wave
	DUPLICATE/O/D radius d_lnr_2, int_Kc_n_n_J_sum, jacobian//, dn_dt_ind
	DUPLICATE/O/D radius dn_dt_ind_raw
	VARIABLE i, m, b, a1, int_Kc_n_n_J, n_avg_a1, d_lnr_1
	DUPLICATE/O/D radius upper, fraction,lower
	a1 = radius[bin]			
	a2 = radius		  			// defines a2 in same order as radius; from large to small
	//a2 = radius[L-1-p]				// defines a2 in reverse order: from small to large size
	// Identical results if order gets reversed: change all commands where it sais " from small to large size".
	a = ( (a1)^3 + (a2)^3 )^(1/3)		// defines a conform (a1)^3 + (a2)^3 = (a)^3
	n_avg_a2 = n_avg				//n_avg_a2 = n_avg[L-1-p] 	// from small to large size
	//n_avg_a2 = (n[p][scan]+n[p][scan+1]) /2 // equivalent to using n_avg, but this way it is 
	n_avg_a1 = n_avg[bin]			// variable // usable by back_calc as well (correct_losses(bin)).
	//n_avg_a1 = (n[bin][scan]+n[bin][scan+1]) /2
	d_lnr_1 = d_lnr[bin]
	d_lnr_2 = d_lnr				// from small to large size
	jacobian = (a/a1)^3
	// Three ways of calculating coag_prod_rate: Cases 1, 2, & 3:
	// Case 1: dNc/dt = Int Int Kc n1 n2 d_lnr1 d_lnr2
	// Case 2: dNc/dt = Int Int Kc n1 n2 (1/d_lnr) d_lnr1 d_lnr2 d_lnr
	// Case 3: dNc/dt = Int Int Kc n1 n2 (r/r1)^3 d_lnr2 d_lnr ("regular" way; gives wrong results!)
	// Case 1 & 2 give the same result; they are mathematically identical. 
	// Case 1 & 2 Jacobian is sustituted: (r/r1)^3 = d_lnr1/d_lnr (by differentiating r^3=(r1)^3+(r2)^3)
	
	//dn_dt_ind_raw = coag_rate_constant(a1,a2)*n_avg_a1*n_avg_a2*d_lnr_2 // case 1
	dn_dt_ind_raw = coag_rate_constant(a1,a2)*n_avg_a1*n_avg_a2*d_lnr_2*d_lnr_1 // case 1* & 2
	// This is dN(a)/dt, the change in particle number for size a. 
	//dn_dt_ind_raw = coag_rate_constant(a1,a2)*n_avg_a1*n_avg_a2*d_lnr_2*jacobian // case 3
	
	// Not added up (integrated) at this point, because dn_dt_ind_raw is a wave of dN/dt (or dn/dt; case 3) for different values of "a". 
	// They are however multiplied by the integrator d_lnr_2, so once they have been distributed over the size bins, the resulting values of dn/dt
	// for each particlar value of "a" can directly be added up to obtain the integrals.
	
	upper = binarysearch(radius, a)	// the bin neighbouring a on the large side.
	
	// radius goes from large to small, thus radius[upper] is the larger one.
	// radius goes from small to large , thus radius[upper] is the smaller one.
	fraction = ( a^3-(radius[upper+1])^3 ) / ( (radius[upper])^3-(radius[upper+1])^3 )  // "fraction" of dn_dt_ind_raw goes to the larger of the two neighbouring bins.
	
	dn_dt_ind = 0
	m =L-1					// m=0 for small to large size and m=L-1 for large to small
	IF(a[m]<radius[0])
		DO
			i = upper[m] 			//  the bin neighbouring a[m] on the large side.
			//// dn_dt_ind[i] = dn_dt_ind[i] + ( fraction[m]*dn_dt_ind_raw[m] ) 				// case1 & 3; larger size bin
			//// dn_dt_ind[i+1] = dn_dt_ind[i+1] + ( (1-fraction[m])*dn_dt_ind_raw[m] ) 			// case 1 & 3; smaller size bin
			// dn_dt_ind[i] = dn_dt_ind[i] + ( fraction[m]*dn_dt_ind_raw[m]/d_lnr[i] ) 			// case 2; larger size bin
			// dn_dt_ind[i+1] = dn_dt_ind[i+1] + ( (1-fraction[m])*dn_dt_ind_raw[m]/d_lnr[i+1] ) 	// case 2; smaller size bin
			
			dn_dt_ind[i] = dn_dt_ind[i] + ( fraction[m]*dn_dt_ind_raw[m]/d_lnr[i] ) 			// case 2; smaller size bin
			dn_dt_ind[i+1] = dn_dt_ind[i+1] + ( (1-fraction[m])*dn_dt_ind_raw[m]/d_lnr[i+1] ) 	// case 2; larger size bin
			
			m = m-1				// m = m+1 for small to large size  amd m=m-1 for large to small size
		WHILE(a[m]<radius[0] )		// from small to large size  (m<L-1)	
	ENDIF
	// At the end of the do-while loop, m = first point # where a>r_max (radius[0]). When resulting size "a" is larger than centre radius of largest bin, both 
	// "fraction" and "upper" get messed up. In that case, dn/dt is 100% appointed to the largest bin, i.e. particles are not allowed to be produced outside the 
	// size range. This way, volume (and thus mass) is (more or less) conserved, and the ratio of total loss/production remains 2, which can be checked.
	DO
		//dn_dt_ind[0] = dn_dt_ind[0] + dn_dt_ind_raw[m] 			// case 1 & 3
		dn_dt_ind[0] = dn_dt_ind[0] + dn_dt_ind_raw[m]/d_lnr[0] 		// case 2
		m = m-1												// m = m+1 for small to large size  and m = m-1 for large to small size
	WHILE( m>=0) 											// ( m<L )	for small to large size and ( m>=0 ) for large to small size
	KILLWAVES a2, a, n_avg_a2, Kc_n_n_J, int_Kc_n_n_J_wave, upper, fraction
	KILLWAVES d_lnr_2, int_Kc_n_n_J_sum, jacobian, dn_dt_ind_raw //, dn_dt_ind
END


FUNCTION coag_rate_constant(R1,R2)
	// kinetic correction is included (enhancement factor). calculates the coagulation rate constant for a pair of particles.
	// cgs units: K_coag in cm^3 /s ,  equations given in CHEM5740-First Handout are used: convention for identical particles.
	// "factor of two" included in rate expression, not in rate constant K_coag.
	VARIABLE R1, R2 
	VARIABLE K_coag, Kc_kin, Kc_con, c_ij, red_mass, a1, a2, D_ij, t_ij, beta_const
	VARIABLE xx, yy, E0, Einf, RedHam
	NVAR Boltz, Temp, density
	a1 = 1.0e-7*R1                        													// convert from nm to cm
	a2 = 1.0e-7*R2  
	red_mass = (((4/3)*density*PI)^2 *(a1^3)*(a2^3)) / (((4/3)*density*pi)*((a1^3)+(a2^3))) 		// reduced mass (g)
	c_ij = SQRT(8*Boltz*Temp/(PI*red_mass))											// mean relative speed of the pair of particles (cm/s)
	// compute the enhancement factors [Chan and Mozurkewich, 2001]. 
	//RedHam = 16					 		// reduced Hamaker Constant //ergs  // 1.6E6 Joule
	RedHam = 6.4E-13						 //ergs  // 6.4E-20 Joule
	xx = ln(1+RedHam)
	yy = SQRT(RedHam)
	E0 = 1 + 0.0757*xx + 0.0015*xx^3
	Einf = 1 + ((yy/SQRT(3))/(1+0.0151*yy)) - 0.186*xx - 0.0163*xx^3
	Kc_kin = (a1+a2)^2 * 0.5*PI*c_ij*Einf 		// Kinetic regime rate constant
	D_ij = diffusioncoeff(R1) + diffusioncoeff(R2)
	Kc_con = 2*PI*D_ij*(a1+a2)*E0 			// Continuum regime rate constant
	t_ij = Kc_kin/(2*Kc_con)       				// same as tij/2 in Sceats
	beta_const = SQRT(1+(t_ij^2)) - t_ij
	K_coag = Kc_kin*beta_const				// Transition regime rate constant (Sceats, 1989)
	
	//IF(R1<R2) 
	//	RETURN 1.5*K_coag * (1-(2/R1))
	//ELSE 
	//	RETURN 1.5*K_coag * (1-(2/R2))
	//ENDIF
	//if (R1<4.6)%| (R2<4.6)
	//	RETURN K_coag * 2
	//ELSEIF (R1<5)%| (R2<5)
	//	RETURN K_coag * 1.5
	//ELSE
	RETURN K_coag    
	//ENDIF
END  
  

FUNCTION first_order_rate()
	NVAR L
	WAVE d_lnr, n_avg
	MAKE/O/D/N=(NUMPNTS(n_avg)) n_dlnr, first_rate
	n_dlnr = n_avg * d_lnr
	first_rate = SUM(n_dlnr, 0, p)
	KILLWAVES n_dlnr
END


FUNCTION mass_particle(R)
	VARIABLE R							// [nm]
	VARIABLE particle_volume
	NVAR density							// density of the particle [g/cc]
	particle_volume = (4/3) * PI * (R*1e-7)^3	// [cc] or [cm^3]
	RETURN particle_volume*density		// [g]
END



	
//FUNCTION Call_Fitting_dNc_dt(ctrlName) : ButtonControl                     
 	//STRING ctrlName
FUNCTION Call_Fitting_dNc_dt()
	Graph_fitting_dNc_dt()
END 


FUNCTION Graph_fitting_dNc_dt() : graph 
	NVAR L
	PauseUpdate; Silent 1		// building window...
	Display /W=(2.25,80,283.5,296) dNc_dt_calc vs radius as "fitting dNc/dt"
	AppendToGraph dNc_dt vs radius
	AppendToGraph/L=Res_Left rel_res vs radius
	AppendToGraph/L=Res_Left_1 Res_dNc_dt vs radius
	ModifyGraph mode(dNc_dt)=3,mode(rel_res)=3,mode(Res_dNc_dt)=3
	ModifyGraph marker(dNc_dt)=41
	ModifyGraph lSize(dNc_dt_calc)=1.5
	ModifyGraph lStyle(dNc_dt_calc)=2,lStyle(Res_dNc_dt)=2
	ModifyGraph rgb(dNc_dt_calc)=(0,0,0),rgb(dNc_dt)=(65280,0,0),rgb(Res_dNc_dt)=(65280,0,0)
	ModifyGraph grid(Res_Left)=2,grid(Res_Left_1)=1
	ModifyGraph log(bottom)=1
	ModifyGraph tick=2
	ModifyGraph mirror=1
	ModifyGraph nticks(Res_Left)=2,nticks(Res_Left_1)=2
	ModifyGraph minor(left)=1,minor(bottom)=1
	ModifyGraph sep(left)=15,sep(bottom)=10,sep(Res_Left_1)=2
	ModifyGraph highTrip(Res_Left)=5
	ModifyGraph standoff=0
	ModifyGraph axOffset(bottom)=-0.8
	ModifyGraph lblPos(left)=49,lblPos(Res_Left)=45,lblPos(Res_Left_1)=34
	ModifyGraph lblLatPos(Res_Left)=2
	ModifyGraph freePos(Res_Left)=0
	ModifyGraph freePos(Res_Left_1)=0
	ModifyGraph axisEnab(left)={0,0.525}
	ModifyGraph axisEnab(Res_Left)={0.81,1}
	ModifyGraph axisEnab(Res_Left_1)={0.5725,0.7625}
	Label left "dNc/dt (cm\\S-3\\Ms\\S-1\\M)"
	Label bottom "radius (nm)"
	Label Res_Left "rel. residuals"
	Label Res_Left_1 "abs. residuals"
	SetAxis/A/N=1 left
	SetAxis/A/N=1 bottom
	SetAxis/A/N=1/E=2 Res_Left
	SetAxis/A/N=1/E=2 Res_Left_1
	//ErrorBars dNc_dt Y,wave=(sd_dNc_dt,sd_dNc_dt)
	Legend/N=text0/J/F=0/B=1/A=MC/X=25.45/Y=-8.74 "\\s(dNc_dt_calc) dNc/dt fitted\r\\s(dNc_dt) dNc/dt measured"
	AppendText "\\s(rel_res) rel_res\r\\s(Res_dNc_dt) Res_dNc_dt"
	Cursor/P A dNc_dt 0 ;Cursor/P B dNc_dt L
	ShowInfo
END


//FUNCTION Call_Fit_range(ctrlName) : ButtonControl       
 	//STRING ctrlName 
FUNCTION Call_Fit_range() 
 	NVAR K
 	VARIABLE first_g,last_g
 	first_g =0
	last_g = K
	Fit_range(first_g,last_g)  	
END     
  
      
FUNCTION Fit_range(first_g, last_g)
	// Performs the fit for a specified number of time intervals. Those time intervals need to be stored in one matrix.
	// Specify value for "number" in beginning and end for the range to be fit. Fitting parameters and chi-squared values written to parameters_matrix.
	// Weighted residuals (res/sd) written to weighted_residuals_matrix. first and last scan (time interval) to be analyzed
	
	VARIABLE first_g, last_g 	// first and last scan (time interval) to be analyzed
	DoWindow/F Graph_fitting_dNc_dt
	VARIABLE i, scan, m 		// i for the row index of w_all (=parameter)
							// scan for the column index (=time interval), m for the row index of res_all (=bin number)
	WAVE w, bin, sd_dNc_dt, dNc_dt, dNc_dt_calc, coag_loss_rate //, Res_dNc_dt, rel_res
	WAVE dn_dt_coag_loss, dn_dt_coag_prod
	DUPLICATE/O/D bin rel_res, res_dNc_dt
	DUPLICATE/D/O w w_sigma
	NVAR L, K 
	DUPLICATE/O/D dNc_dt UC_dNc_dt, LC_dNc_dt, UP_dNc_dt, LP_dNc_dt
	MAKE/D/O/N=((NUMPNTS(w)+1),(K-1)) w_all
	MAKE/D/O/N=((NUMPNTS(w)),(K-1)) sigma_w
	MAKE/D/O/N=(L,(K-1)) res_all, coag_net_rate //, abs_res_all
	//MAKE/D/O/N=(L,(K-1)) dNc_dt_all 
	MAKE/D/O/N=(K-1) growth, chisq, coag_total_loss_rate, sigma_g, rate_c, Conf_IV_g
	MAKE/T/O/N=1 T_Constraints
	
	//T_Constraints = {"K3>0.1"}//,"K6>0.1"} 					// when coag and wall loss are fitted simultaneously
	//T_Constraints = {"K1>0.7","K1<1"}						// when alpha is fitted
	//T_Constraints = {"K1>0.7","K1<1","K2>0.08","K2<0.12"} 	// when DiffCoeff is fitted
	//T_Constraints = {"K8>0.0001"}//,"K9>20","K9<40"} 		// when svp & surface tension are fitted
	w_all = NaN
	w_all[0][] = -999										// set growth rate to -999 instead of NaN, to not interfere with interpolation
	sigma_w = NaN
	res_all = NaN
	coag_total_loss_rate = NaN
	//coag_total_loss_r_ex_nano = NaN
	coag_net_rate = NaN
	//dNc_dt_all = NaN
	res_dNc_dt = NaN
	Conf_IV_g = NaN
	scan = first_g
	DO
		WAVE W_ParamConfidenceInterval	
		start (scan)
		losses( )		
		UC_dNc_dt = NaN 
		LC_dNc_dt = NaN 
		UP_dNc_dt = NaN 
		LP_dNc_dt = NaN
		dNc_dt_calc =NaN 
		
		coag_total_loss_rate[scan] = coag_loss_rate[L-1]					// dN(total)/dt due to coag loss
		//coag_total_loss_r_ex_nano[scan] = coag_loss_rate[29] 			// idem, excl nanoDMA (Pacific 2001)
		coag_net_rate[][scan] = dn_dt_coag_loss[p] - dn_dt_coag_prod[p] 	// dn/dt due to coag loss-prod, for each size bin.
		   
		// performs the fit, to be repeated for each time interval specified.
		//FuncFit /Q/H="01111111" GDE w dNc_dt[first_g, last_g] /I=1 /W=sd_dNc_dt /D=dNc_dt_calc /R ///C=T_Constraints
		//FuncFit /Q/G/H="01111111" GDE w dNc_dt[first_g, last_g] /W=sd_dNc_dt /I=1 /D=dNc_dt_calc /R /F={0.95, 7, Contour, UC_dNc_dt, LC_dNc_dt, UP_dNc_dt, LP_dNc_dt}
	        //FUNCFIT/Q/G/H="01111111" GDE w dNc_dt[first_g, last_g] /W=sd_dNc_dt /I=1 /D=dNc_dt_calc /R /F={0.95, 4} //7, Contour, UC_dNc_dt, LC_dNc_dt, UP_dNc_dt, LP_dNc_dt}
		// FuncFit /Q/H="01111111" GDE w dNc_dt /I=1/W=sd_dNc_dt /D=dNc_dt_calc /R //C=T_Constraints

		//FUNCFIT/Q/G/H="01111111" GDE w dNc_dt[first_g, last_g] /W=sd_dNc_dt /D=dNc_dt_calc /R /F={0.95, 4}
		FuncFit /Q/G/H="01111111" GDE w dNc_dt[first_g, last_g] /W=sd_dNc_dt /I=1 /D=dNc_dt_calc /R /F={0.99, 4}//, Contour, UC_dNc_dt, LC_dNc_dt, UP_dNc_dt, LP_dNc_dt}
			
		rel_res = Res_dNc_dt / sd_dNc_dt 	   					// chi = relative residual / standard deviation
		 
		i = 0  
		DO    
			w_all[i][scan] = w[i]   			 					// save fit parameters in matrix w_all
			sigma_w[i][scan] = w_sigma[i]	 					// save estimated error of fit parameters in matrix sigma_w
			i = i + 1 
			m = 0  
			DO 
				res_all[m][scan] = Res_dNc_dt[m] / sd_dNc_dt[m] 	// matrix of rel_res // during the fit, weighted residuals (res./sd.) are written to matrix res_all
				// abs_res_all[m][scan] = Res_dNc_dt[m]
				// dNc_dt_all[m][scan] = dNc_dt[m]				// during the fit, the values for dNc_dt are written to matrix dNc_dt_all				
				m = m+1
			WHILE(m < L)
		WHILE(i < NUMPNTS(w))
		w_all[NUMPNTS(w)][scan] = V_chisq	
		// save chi-squared of fit
		Conf_IV_g[scan] = W_ParamConfidenceInterval[0]		// confidence interval of fit parameter w[0] = growth rate
		// sigma_w[numpnts(w)+1][scan] = w_sigma[0]				// save st.dev. of growth estimate
		scan = scan + 1
	WHILE(scan <= last_g)		 								// defines the number of the last time interval to be analyzed
	growth = w_all[0][p]	  									// wave of size-averaged growth rate (at infinite Kn) for each time interval
	i = 0
	DO
		IF(growth[i] == -999 || growth > 300 || growth < -300 )
			growth[i] = NaN 
		ENDIF
		i = i+1 
	WHILE(i < NUMPNTS(growth))   
	
	rate_c = w_all[8][p] * 3600									// wave of size-averaged rate (at infinite Kn) for each time interval
	//duplicate/O/D growth g_smooth
	//Smooth 7, g_smooth // adapt smoothing according to amount of noise in growth vs. time
	
	chisq = w_all[NUMPNTS(w)][p]									// wave of chi-squared values
	sigma_g = sigma_w[0][p]									// wave of IGOR reported (estimated) sd of growth estimate
	duplicate/O/D sigma_g error_g_all, error_scatter 				// V_npnts is number of size bins included in the fit
	error_scatter = SQRT(chisq/(V_npnts-1)) 						// error related to scatter around the obtained value for fit parameter
	
	// This assumes only one parameter being fitted (# degrees of freedom = V_npnts - # fit parameters)
	//error_g_all = sigma_g * error_scatter 						// combined error (estimated (expected) and from scatter)
	error_g_all = Conf_IV_g * error_scatter	 					 // combined error (estimated (expected) and from scatter)
	//killwaves error_scatter, sigma_g, w_sigma
	//separate_factors()
	//chi_squared(first, last) // Elaborate error analysis using chi-values	
END
   


 
FUNCTION GDE(w,bin) : FitFunc		// describes fitting function: General Dynamic Equation in cumulative form, n=dN/dln(r)
	WAVE w						// wave holding the parameters that can be fitted
	VARIABLE bin
	NVAR L
	WAVE Nc_avg, n_avg, radius, coag_loss_rate, nn_avg2, nn_avg
	WAVE wall_diff_rate, wall_grav_rate, coag_prod_rate	// , first_rate,appear,appear_rate	// Input data tables.
	VARIABLE growth, wall, dNc_dt_calc, coag, rate		//two_dim_avg(w)
	// VARIABLE appear	 		// appearance rate of particles at smallest size (due to nucleation)
	
	// growth applies Fuchs/Sutugin correction for all sizes.
	// w[0] is the radius growth rate (nm/h) divided by gamma at infinite Knudsen number (g_0)
	// w[1] is alpha
	// w[2] is the diffusion coefficient of the condensing species 
	// w[3] is coagulation loss and production multiplier, which multiplies both by the same factor
	// w[4] is coagulation loss multiplier
	// w[5] is coagulation production multiplier
	// w[6] is diffusional wall loss multiplier
	// w[7] is gravitational wall loss multiplier
	// w[3] should not be fitted simultaneously with either w[4] or w[5]
	// the "multiplier" parameters are set equal to 1 to give the calculated value
	// growth = ((w[0]-evaporation(radius[bin], w[8], w[9], 1))/3600)*nn_avg[bin]*Gamma_meas(radius[bin], w[1], w[2]) 
	
	growth = (w[0]/3600)*nn_avg[bin]*Gamma_meas(radius[bin], w[1], w[2]) 
	coag = w[3] * ( (coag_loss_rate[bin]*w[4])  - (coag_prod_rate[bin]*w[5]))
	wall = (wall_diff_rate[bin]*w[6]) + (wall_grav_rate[bin]*w[7]) 	// + w[8]*Nc_avg[bin] 
	// rate = first_rate[bin] * w[8]
	// appear = appear_rate[bin] * w[9]
	//IF(bin == (L-1) )
	//	appear = w[9]
	//ELSEIF(bin < (L-1) )
	//	appear = 0
	//ENDIF

	dNc_dt_calc = growth - coag - wall  //+ rate// + appear
	RETURN dNc_dt_calc 
END

FUNCTION Gamma_meas( radius, alpha, DiffCoef)   // Returns the apparent sticking coefficient for the given parameters.
	// DiffCoef is the gas phase diffusion coeffcient of the condensing species  ==  0.1 cm^2/s (+/- 10%) for H2SO4
	VARIABLE radius, alpha, DiffCoef
	VARIABLE lambda, resistance
	NVAR speed_gas
	lambda = (1e7*3*DiffCoef) / speed_gas		// lambda (nm) of condensible species
	resistance = (1/alpha) + ((3*radius)/(4*lambda)) - ((0.47*radius)/(radius+lambda))
	RETURN 1/resistance 						// dimensionless gamma
END


FUNCTION Call_Graph_growth_vs_time(ctrlName) : ButtonControl     
 	STRING ctrlName
       	 Graph_growth_vs_time( )
END

FUNCTION save_GR_Jnuc(ctrlName) : ButtonControl  
       STRING ctrlName
       WAVE t_midgrowth,growth_smth,t_meas_long,Jnuc
       SMOOTH/B=2 5, J_nucl
	Save/J/M="\r\n" t_midgrowth, growth_smth as "GR_.txt"   
	Save/J/M="\r\n" t_form_long,J_nucl as "Jnucl_.txt"   
END	
  
FUNCTION Graph_growth_vs_time() : Graph
	//Save/J/M="\r\n" t_midgrowth as "GR_time.txt"   
	//Save/J/M="\r\n" growth as "GR.txt"    
	wave growth
	IF(growth > 300 || growth < -300)
		growth=NAN    
	ENDIF  
	PauseUpdate; Silent 1		 
	Display /W=(294.75,333.5,590.25,552.5) growth vs t_midgrowth as "growth_vs_time"
	ModifyGraph mode=3
	ModifyGraph marker=19
	ModifyGraph tick=2    
	ModifyGraph mirror=1 
	ModifyGraph minor=1
	ModifyGraph sep(left)=15,sep(bottom)=10
	ModifyGraph lblMargin(bottom)=2
	ModifyGraph standoff=0
	ModifyGraph axOffset(left)=-2.71429,axOffset(bottom)=-0.1875
	ModifyGraph lblLatPos(left)=-1,lblLatPos(bottom)=-9
	ModifyGraph dateInfo(bottom)={1,1,0}
	Label left "growth rate (nm/h)"
	Label bottom "time"
	SetAxis/A/N=1 left 
	ModifyGraph dateInfo(bottom)={1,0,0}
	
	print "GR=", MedianEO(growth, 0, DIMSIZE(t_midgrowth,0)-1)
	print "GR=", Median(growth, 0, DIMSIZE(t_midgrowth,0)-1)
	
END
  
//###################################################################
//################### Obtaining Nucleation Rates ##########################
//###################################################################

FUNCTION Nucl_Rate(ctrlName) : ButtonControl      
 	STRING ctrlName
 	NVAR K
 	VARIABLE first_g,last_g
 	first_g = 0
	last_g = K
	Fit_Kw() 	
	Fit_kc1()  
	Fit_Kc()
	back_calc_mult(first_g, last_g)
END
  
 	 
FUNCTION Fit_Kw()   // fit the wall loss rate constant to powerlaw of form K_w = wall_const/radius
	MAKE/D/O/N=20 a_small, Kw_s, Kw_fit_s
	VARIABLE/G wall_const
	VARIABLE i
	WAVE w      
	NVAR r_nucl 
	a_small[0] = r_nucl
	i = 1
	DO
		a_small[i] = a_small[i-1] + 0.5
		i = i+1
	WHILE(i<20)
	Kw_s = wall_loss_diff(a_small) * w[6]
	K0 = 0
	K1 = 0.00098*w[6]
	K2 = -1
	//CurveFit/G Power Kw_s /X=a_small /D=Kw_fit_s 
	CurveFit/G/H="101" Power Kw_s /X=a_small /D=Kw_fit_s 
	wall_const = K1					// nm/s
	//KILLWAVES/Z a_small, Kw_s, Kw_fit_s
END


FUNCTION Fit_kc1() 	
	// Fitting a power-law to the pseudo first order coagulation rate constant to be used for back-calculation. 
	// Kc1(a) = Sum(Kc2(a,A)*n(A))	 = i+C*a^p, where a<A.  Fitting parameters are f(t) (a function of time 
	// ("scan"), because n(A) is.).  Scavenging by larger particles only is indicated by "_sc" extension of wave name.

	WAVE radius, d_lnr, dNdlogDp, w
	DUPLICATE/O/D radius Kc2_r_var, Kc1_wave, Kc1_wave_sc, Kc1_wave_sc_acc
	NVAR L, K, r_nucl
	VARIABLE n_a_small
	n_a_small = 20
	MAKE/D/O/N=(n_a_small) a_small, Kc1, Kc1_fit, Kc1_sc, Kc1_fit_sc, Kc1_sc_acc, Kc1_fit_sc_acc
	MAKE/D/O/N=(K) i_kc1_fit, c_kc1_fit, p_kc1_fit, kc1_2nm, kc1_2nm_fit
	MAKE/D/O/N=(K)  i_kc1_fit_sc, c_kc1_fit_sc, p_kc1_fit_sc, kc1_2nm_sc, kc1_2nm_fit_sc, kc1_2nm_sc_acc, kc1_2nm_fit_sc_acc
	MAKE/D/O/N=(K)  i_kc1_fit_sc_acc, c_kc1_fit_sc_acc, p_kc1_fit_sc_acc
	MAKE/D/O/N=(K)  kc1_8nm, kc1_8nm_fit, kc1_8nm_sc, kc1_8nm_fit_sc, kc1_8nm_sc_acc, kc1_8nm_fit_sc_acc
	MAKE/D/O/N=(3,K) w_Kc1_fit, w_Kc1_fit_sc, w_Kc1_fit_sc_acc
	MAKE/D/O/N=3 w_Kc1_fit1, w_Kc1_fit1_sc, w_Kc1_fit1_sc_acc
	VARIABLE i, r_var, scan, j
	MAKE/O/T/N=4 T_Constraints_coag
	a_small[0] = r_nucl
	i = 1 
	DO
		a_small[i] = a_small[i-1] + 0.5
		i = i+1
	WHILE(i<n_a_small)
	T_Constraints_coag[0] = {"K0>1e-15", "K1>1e-15", "K2<-1", "K2>-2"}
	scan = 0 				// the parameter for scan number, i.e. time
	DO
		i = 0				 // the parameter for a_small
		DO
			Kc1_wave = 0 // not NaN, because then Kc1_wave[L-1] returns NaN.
			Kc1_wave_sc = 0
			r_var = a_small[i]
			Kc2_r_var = w[3] * coag_rate_constant(r_var,radius) // 2nd order K_coag
			Kc1_wave = dNdlogDp[p][scan] * Kc2_r_var * d_lnr	// ps. 1st order K_coag for sc. by all part.
			j = 0			 // the parameter for bin number of scavenging particle
			DO 			 // only inlcude scavenging by larger particles than r_var itself:
				Kc1_wave_sc[j] = dNdlogDp[j][scan] * Kc2_r_var[j] * d_lnr[j] // ps. 1st order K_coag
				j = j+1
			WHILE( (j<L) && (radius[j]>r_var) )
			j = 0 		// the parameter for bin number of scavenging particle
			DO			// only inlcude scavenging by much larger, accumulation mode particles:
				Kc1_wave_sc_acc[j] = dNdlogDp[j][scan] * Kc2_r_var[j] * d_lnr[j] // ps. 1st order K_coag
				j = j+1
			WHILE( (j<21) && (radius[j]>r_var) )
			//Kc1[i] = sum(Kc1_wave, 0, L) 		// identical to rectangular integration
			INTEGRATE Kc1_wave 			// summation over all bins: rectangular integration
			INTEGRATE Kc1_wave_sc			// summation over all bins: rectangular integration
			INTEGRATE Kc1_wave_sc_acc 		// summation over all bins: rectangular integration
			Kc1[i] = Kc1_wave[L-1]
			Kc1_sc[i] = Kc1_wave_sc[L-1]
			Kc1_sc_acc[i] = Kc1_wave_sc_acc[L-1] // note that Kc1 has a different meaning here than in previous function.
			i = i+1
		WHILE(i<n_a_small) 
		w_Kc1_fit1_sc[0] = 1e-7
		w_Kc1_fit1_sc[1] = 1e-5
		w_Kc1_fit1_sc[2] = -1.7
		CurveFit/G Power KwCWave=w_Kc1_fit1_sc Kc1_sc /X=a_small /D=Kc1_fit_sc /C=T_Constraints_coag 
		w_Kc1_fit_sc[][scan] = w_Kc1_fit1_sc[p]
		Kc1_2nm_sc[scan] = Kc1_sc[3] 		// to compare fitted with "real" ps first order Kc. (coag with all part.)
		Kc1_8nm_sc[scan] = Kc1_sc[15] 	// to compare fitted with "real" ps first order Kc. (coag with all part.)
		w_Kc1_fit1[0] = 1e-7
		w_Kc1_fit1[1] = 1e-5
		w_Kc1_fit1[2] = -1.7
		CurveFit/G Power KwCWave=w_Kc1_fit1 Kc1 /X=a_small /D=Kc1_fit /C=T_Constraints_coag 
		w_Kc1_fit[][scan] = w_Kc1_fit1[p]
		Kc1_2nm[scan] = Kc1[3] 			// to compare fitted with "real" ps first order Kc. (sc. by larger part. only)
		Kc1_8nm[scan] = Kc1[15]			// to compare fitted with "real" ps first order Kc. (sc. by larger part. only)
		
		//CurveFit/G Power Kc1 /X=a_small1 /D=Kc1_fit
		//w_Kc1_fit1_sc_acc[0] = 1e-7
		//w_Kc1_fit1_sc_acc[1] = 1e-5
		//w_Kc1_fit1_sc_acc[2] = -1.7
		//CurveFit/G Power KwCWave=w_Kc1_fit1_sc_acc Kc1_sc_acc /X=a_small /D=Kc1_fit_sc_acc /C=T_Constraints_coag 
		//w_Kc1_fit_sc_acc[][scan] = w_Kc1_fit1_sc_acc[p]
		//Kc1_2nm_sc_acc[scan] = Kc1_sc_acc[3] // to compare fitted with "real" ps first order Kc. (coag with all part.)
		//Kc1_8nm_sc_acc[scan] = Kc1_sc_acc[15] // to compare fitted with "real" ps first order Kc. (coag with all part.)
		scan = scan+1
	WHILE(scan<K) // 25)
	i_kc1_fit_sc = w_Kc1_fit_sc[0][p]
	c_kc1_fit_sc = w_Kc1_fit_sc[1][p] 
	p_kc1_fit_sc = w_Kc1_fit_sc[2][p]
	//i_kc1_fit_sc_acc = w_Kc1_fit_sc_acc[0][p]
	//c_kc1_fit_sc_acc = w_Kc1_fit_sc_acc[1][p]
	//p_kc1_fit_sc_acc = w_Kc1_fit_sc_acc[2][p]
	i_kc1_fit = w_Kc1_fit[0][p]
	c_kc1_fit = w_Kc1_fit[1][p]
	p_kc1_fit = w_Kc1_fit[2][p]
	kc1_2nm_fit = i_kc1_fit + c_kc1_fit * 2^p_kc1_fit
	kc1_2nm_fit_sc = i_kc1_fit_sc + c_kc1_fit_sc * 2^p_kc1_fit_sc
	//kc1_2nm_fit_sc_acc = i_kc1_fit_sc_acc + c_kc1_fit_sc_acc * 2^p_kc1_fit_sc_acc
	kc1_8nm_fit = i_kc1_fit + c_kc1_fit * 8^p_kc1_fit
	kc1_8nm_fit_sc = i_kc1_fit_sc + c_kc1_fit_sc * 8^p_kc1_fit_sc
	//kc1_8nm_fit_sc_acc = i_kc1_fit_sc_acc + c_kc1_fit_sc_acc * 8^p_kc1_fit_sc_acc
	//KILLWAVES/Z Kc2_r_var, Kc1_wave, a_small, Kc1, Kc1_fit, w_Kc1_fit, w_Kc1_fit1, T_Constraints_coag
END 



FUNCTION Fit_Kc()
	// Fitting a power-law to the coagulation rate constant of 2 particles of different size, to be used for back-calculation. 
	// Kc(a,A) = i+c*a^p	where a<A, and fitting parameters are f(A). Creates a size dependent wave for each fitting 
	// parameter: intercept (i_fit_c); constant (c_fit_c); power (p_fit_c)
	WAVE radius, w
	VARIABLE i
	NVAR L, r_nucl
	MAKE/D/O/N=20 a_small, Kc_x, Kc_fit_x
	MAKE/D/O/N=3 w_Kc_fit1
	MAKE/D/O/N=(3,L) w_Kc_fit
	DUPLICATE/D/O radius c_fit_c, i_fit_c, p_fit_c, Kc_5, Kc_5_fit
	a_small[0] = r_nucl
	i = 1
	DO
		a_small[i] = a_small[i-1] + 0.5
		i = i+1
	WHILE(i<20) 			// a_small from 0.5 to 10 nm
	w_Kc_fit[0][] = 0
	w_Kc_fit[1][] = 1e-6
	w_Kc_fit[2][] = -1.7
	K0 = 0;K1 = 1e-6;K2 = -1.8;
	Make/O/T/N=1 T_Constraints_coag
	T_Constraints_coag[0] = {"K0 > 0"}
	i = 0
	DO 					// Calculate K_coag as it is actually used in the fit procedure by including w[3]
		Kc_x = coag_rate_constant(a_small, radius[i]) * w[3]
		w_Kc_fit1 = w_Kc_fit[p][i]
		CurveFit/G Power KwCWave=w_Kc_fit1 Kc_x /X=a_small /D=Kc_fit_x /C=T_Constraints_coag 
		w_Kc_fit[][i] = w_Kc_fit1[p]
		//w_Kc_fit contains fitting parameters for Kc as a function of radius(A).
		i = i+1
	//WHILE(radius[i] > a_small[19])
	WHILE(i<L)
	i_fit_c = w_Kc_fit[0][p]
	c_fit_c = w_Kc_fit[1][p]
	p_fit_c = w_Kc_fit[2][p]
	Kc_5 = w[3]*coag_rate_constant(5, radius)
	Kc_5_fit = i_fit_c + c_fit_c * 5^p_fit_c
	//KILLWAVES Kc_x, Kc_fit_x, w_Kc_fit, w_Kc_fit1, T_Constraints_coag //, a_small
END



FUNCTION back_calc_mult(first_g, last_g)  
	// Creates matrix of radii for series of scans (measurement times) with radii calculated backwards in time due to 
	// (smoothed) growth rate. Creates matrix of formation times (row denotes bin,  column denotes measurement time ("scan")
	VARIABLE first_g, last_g
	VARIABLE bin
	NVAR r_nucl, L, K		 // number 
	MAKE/O/D/N=(K,L) r_back
	MAKE/O/D/N=(L,K) t_form, r_meas_long, nn_form_cor
	VARIABLE t, t_form_ind_var, r_var_init, i, j, m, scan
	WAVE growth, radius, w, t_wave, t_end_wave, w, nn, dNdlogDp,dt
	//DUPLICATE/O/D t_wave r_back_ind //, t_midnucl
	MAKE/O/D/N=(NUMPNTS(t_wave)) r_back_ind
	//DUPLICATE/O/D growth delta_r_back, dt_wave //, t_avg_wave, t_midgrowth//, growth_orig
	MAKE/O/D/N=(NUMPNTS(growth)) delta_r_back, dt_wave //, t_avg_wave, t_midgrowth//, growth_orig
	
	//Smooth 3, growth // adapt smoothing according to amount of noise in growth vs. time. NOTE that smoothing is useful for 
	// noisy (atmospheric) data, but often not desirable for smog-chamber data. 
	DUPLICATE/O/D radius radius_back, delta_t_back, nn_form_cor_ind //, t_form_ind
	MAKE/D/O/N=(L,K) dep_F, coagsc_F //, wmcoag_F, b1_F, wmcpulse_F, b1p_F
	//make/O/D/N=(K*L,K) n_back_all, r_back_all 		// only needed when reconstructing size distributions
	MAKE/O/D/N=(K*L) n_back_t1, r_back_t1
	WAVE t_midnucl, t_midgrowth

	dt_wave = t_wave[p+1] - t_wave[p] 			// wave of time intervals (seconds)
	//dt_wave = dt
	  
	//t_avg_wave = (t_wave[p]+t_wave[p+1]) / 2
	//t_midnucl = t_wave + 40 // VALID FOR 11-12 ONLY ! ! !  and 11-19. and 11-16.
	//i = 0
	//do	// ADAPT FOR EACH CASE STUDY.
	//	t_midnucl[i] = t_wave[i] + 115 
	//	i = i+1
	//while (i<=55)
	//t_midgrowth = (t_midnucl[p]+t_midnucl[p+1]) / 2 // avg time between the maxima of 2 subsequent dist's.
	// Use t_midnucl for "point-data" such as J, n, N. Use t_midgrowth for "time interval-data" such as G. 
	
	r_back= 999  // instead of NaN, because interpolation can't handle NaN.
	// NOTE that use of "999" will interfere if cubic or smooth spline interpolation is used.
	
	t_form = 999
	dep_F = NaN
	coagsc_F =NaN
	//wmcoag_F = NaN
	//wmcpulse_F = NaN
	//b1_F = NaN
	//n_back_all = NaN
	//r_back_all = NaN
	//t_form_ind = NaN
	nn_form_cor = NaN
	nn_form_cor_ind = NaN
	delta_r_back = growth*dt_wave/3600  	// delta_r_back[scan] valid for change from n[scan] to n[scan+1]
	// these waves are f(t) at measurement time defined by "scan" (minimum value 1).
	//t = first_g // t = timestep of backwards calculation; starts from "scan", to be counted down.
	
	scan = first_g 					// scan = number of scan to be back-calculated from
	DO 
		//t = scan 
		//r_back[scan][] = radius[p]
		//bin = 0					// starting at 406 nm radius
		bin = 0 //90 //44 //7		// starting at smaller radius
		//print "Initial bin =",bin		
		DO
			t = scan			
			r_back[scan][bin] = radius[bin]
			//r_var_init = radius[bin]
			DO	// One column of the matrix r_back gets calculated...
				r_back[t-1][bin] = r_back[t][bin] - (delta_r_back[t-1]*gamma_meas(r_back[t][bin], w[1], w[2]))
				// matrix of backwards calculated evolution of radii for one specific scan number.
				t=t-1			
			 WHILE(t>0) 
			 
			// this column gets written to the wave r_back_ind.
			r_back_ind = r_back[p][bin] // wave of backwards radius for one particular
			// value of bin (gets overwritten for each new looping bin-value) in order to
			// interpolate where the radius passes the critical value denoted by r_nucl
				
			//t_form_ind_var = INTERP(r_nucl, r_back_ind, t_midnucl) 		// formation time, for nucleated particles around 13 nm.
			t_form_ind_var = INTERP(r_nucl, r_back_ind, t_wave) 		// formation time, if nucleated part are small.
			//t_form_ind_var = INTERP(r_nucl, r_back_ind, t_avg_wave)	// if nucleated part are large.
			//if ( t_form_ind_var < t_wave[0] ) 
			 
										 
			IF( t_form_ind_var < t_wave[0] )   
				t_form[bin][scan] = 999 			// discard formation times that are earlier than measurements began.
			//ELSEIF( t_form_ind_var > t_wave[K-1] ) 	// This "can not" happen, but strangely enough, it does.
				//t_form[bin][scan] = 999 			// discard formation times that are later than measurements ended.
			ELSEIF( t_form_ind_var > t_wave[scan] ) 	// This "can not" happen, but strangely enough, it does.
				t_form[bin][scan] = 999		// discard formation times that are later than actual measurement
			ELSE   
				t_form[bin][scan] = t_form_ind_var 	// wave of formation times for particular meas. time (scan)				
				//For each value of "bin", the wave "n_back_ind" is overwritten.
				nn_form_cor[bin][scan] = correct_losses_mult(bin,scan)	 // corrected for losses.	
				//fill_n_back(bin, scan)
				//Only executed t_form is later than first measurement time, unless 
				//you are interested in checking n_back_all
				//nn_form_cor_ind[bin] = n[bin][scan]/r_back_ind[scan] 	// n[bin][scan]/radius[bin] not corrected for losses.
				//n_back_ind is now specified.
			ENDIF  
			 
			//NOTE: enable the next two command if you want to look at the backcalculated waves
			//If not, better put previous command within last "else-endif" loop to save computing time.
			//specifies r_back_all & n_back_all for one value of bin & scan.
		 	//nn_form_cor[bin][scan] = correct_losses_mult(bin, scan) 	// corrected for losses.
			//fill_n_back(bin, scan)			
			bin=bin+1
			//print "Updated bin and L=",bin,L 
		//WHILE(bin<36) 
		WHILE(bin<L)
		scan = scan+1
	WHILE(scan<=last_g) 
	//n_t1 = n[p][12] //[scan] // n(r)_back_t1 gives all backwards calculated distributions,
	//n_back_t1 = n_back_all[p][12] 		// present at the physical time (t1="scan").
	//r_back_t1 = r_back_all[p][12] 		// They should resemble the present distribution closely. 
	// update long table with formation times, nucleation rates, etc:
	
	DUPLICATE/O/D t_form t_form_long
	REDIMENSION/D/N=(K*L) t_form_long
	DUPLICATE/O/D nn_form_cor nn_form_long
	REDIMENSION/D/N=(K*L) nn_form_long 	
	MAKE/O/D/N=(K*L) g_long, J_nucl, t_diff_long, t_meas_long	//, J_new_wmcp, J_nucl_new
	
	//duplicate/O/D nn nn_long
	//redimension/N=(K*L) nn_long	
	//duplicate/O/D wmcoag_F wmc_F_long
	//redimension/N=(K*L) wmc_F_long
	//duplicate/O/D wmcpulse_F wmcp_F_long 
	//redimension/N=(K*L) wmcp_F_long
	//duplicate/O/D b1_F b1_F_long
	//redimension/N=(K*L) b1_F_long  
	//duplicate/O/D b1p_F b1p_F_long
	//redimension/N=(K*L) b1p_F_long
	//duplicate/O/D coagsc_F coagsc_F_long
	//redimension/N=(K*L) coagsc_F_long
	//duplicate/O/D dep_F dep_F_long
	//redimension/N=(K*L) dep_F_long
	//duplicate/O/D b1_F_long ratio_b1_long, coag_total, ratio_b1p_long, dep_coag_F_long
	//dep_coag_F_long = dep_F_long*coagsc_F_long
	 
	g_long = INTERP(t_form_long, t_midgrowth, growth) //growth)		          
	//g_long = interp(t_form_long, t_end_wave, growth) //growth)
	 //interpolated growth rate at formation  time, using end times as x-value for g
  
	J_nucl = nn_form_long*g_long/3600 // Nucleation Rate, NOT corrected
	 
	 FOR(i=0;i<(K*L);i+=1)
	 	IF(J_nucl[i] >5E2 || J_nucl[i] < 1E-4)
	 		J_nucl[i]=NAN
	 	ENDIF
	 ENDFOR    
	//DUPLICATE/O/D g_long growth_smth 
	//SMOOTH/B=3 3, growth_smth  
	//DUPLICATE/O/D J_nucl J_nucl_smth  
	//SMOOTH/B=3 3, J_nucl_smth   
	  
	//J_nucl_new = nn_form_long*g_long/3600 // Nucleation Rate, NOT corrected
	//J_new_wmcp = J_nucl_new*wmcp_F_long/wmc_F_long
	//ratio_b1_long = b1_F_long / (dep_F_long*coagsc_F_long)
	//ratio_b1p_long = b1p_F_long / (dep_F_long*coagsc_F_long)
	//coag_total = wmc_F_long*coagsc_F_long
	//MAKE/D/O ratio_b1_h  
	//HISTOGRAM/B={0, 0.1, 20} ratio_b1_long, ratio_b1_h 
	//DUPLICATE/O/D ratio_b1_h ratio_b1
	//ratio_b1[0] = 0.05  
	//i=1  
	//DO   
	//	ratio_b1[i] = ratio_b1[i-1]+0.1
	//	i=i+1
	//WHILE(i<20)	
	// for losses having occurred before time of measurement. Include measured radius and time of measurement in table for clarity:
	r_meas_long = radius[p]
	REDIMENSION/D/N=(K*L) r_meas_long
	i=0
	j=0
	t_meas_long = NaN
	DO				// this nested do-loop writes a t_meas at each r_meas.
		DO
			t_meas_long[i] = t_wave[j]
			i=i+1
		WHILE( i<((j+1)*L) )
		j=j+1  
	WHILE(j<K) 
	t_diff_long = t_meas_long - t_form_long
	//KILLWAVES r_back_all, n_back_all
	//KILLWAVES dep_F, coagsc_F //, wmcoag_F, b1_F, wmcpulse_F, r_meas_long,b1p_F
	//KILLWAVES b1p_F_long,coagsc_F_long,dep_F_long,dep_coag_F_long,ratio_b1p_long
	
	//KILLWAVES r_back, t_form, nn_form_cor, nn_form_cor_ind // , r_back_ind
	//KILLWAVES delta_r_back, dt_wave, radius_back, delta_t_back//, DF_wave,CF_wave
END
 


FUNCTION correct_losses_mult(bin, scan)
	// Corrects n for losses having occurred between t_form and t_meas.
	// Returns nn at t_form, corrected for those losses, for a particular value of radius, as defined by bin.
	VARIABLE bin, scan
	NVAR r_nucl, K, L, Boltz, Temp, Kw_kin_const, B_kin_const, wall_const, inc
	VARIABLE r_meas, r_var, r_var_1, r_var_2, g_r_var, g_r_nucl, CF, DF
	VARIABLE i, j, m, int_Kc_n_cor, n_back_ind_avg, Kc1_r_var, nn_meas
	WAVE r_back, r_back_ind, radius, dt_wave, d_lnr, growth, w, N_total
	WAVE i_fit_c, c_fit_c, p_fit_c, dNdlogDp, d_lnr, dep_F, coagsc_F //, wmcoag_F, b1_F, wmcpulse_F, b1p_F
	WAVE i_kc1_fit, c_kc1_fit, p_kc1_fit
	DUPLICATE/O/D radius Kc_n_cor, int_Kc_n_cor_wave, int_Kc_n_cor_sum
	DUPLICATE/O/D radius Kc1_int1, Kc1_int2, Kc1_int, K_w_diff_fit, Kc2_r_var, Kc1
	DUPLICATE/O/D growth Kw_cor, Kc_cor, CF_wave, DF_wave
	DUPLICATE/O/D r_back_ind n_back_ind //, b, f_1, f_2, f
	DUPLICATE/O/D radius inner_int, inner_int_1, inner_int_2
	WAVE n_t1, n_t2
	n_t1 = dNdlogDp[p][scan]				// to read in from dN/dlnr matrix (produced by awk)
	n_t2 = dNdlogDp[p][scan+inc]			// dN/dlnr (units of r don't matter)
	replace_NaN_by_num(n_t1, 0)
	replace_NaN_by_num(n_t2, 0)
	//SMOOTH/B=2/F/E=3 7, n_t1	// NOTE THE SMOOTHING!!!   
	//SMOOTH/B=2/F/E=3 7, n_t2

	SMOOTH/B=2 5, n_t1			// NOTE THE SMOOTHING!!!   
	SMOOTH/B=2 5, n_t2

	CF_wave = 1 //NaN
	DF_wave = 1 //NaN
	n_back_ind = NaN    
	//n_back_ind[scan] = n[bin][scan]
	//replace_NaN_by_num(n_back_ind, 0)
	n_back_ind[scan] = n_t1[bin]
	nn_meas = n_back_ind[scan]/ radius[bin] 		//r_back_ind[scan] // radius[bin]
	// remember that r_back_ind[scan] = radius[bin]
	r_meas = radius[bin] //r_back_ind[scan]	 	// the radius from which to start the backwards calculation.
	//NVAR wall_const is resulting value from fit to Kw = wall_const/a  (0.00098)
	// wall_const = Kw_kin_const*sqrt(B_kin_const*Boltz*Temp)
	// wall_const calculated algebraically from Kw equation in limit of a->0
	// Kw = (1/R)*wall_const; sqrt(...) = sqrt(D)
	i = scan-1
	DO		// calculates loss processes integrated backwards to t_form
		IF( r_back_ind[i]>r_nucl )
			r_var_2 = r_back_ind[i+1]			// the later (larger) radius
			r_var_1 = r_back_ind[i]				// the earlier (smaller) radius
			r_var = (r_var_1+r_var_2) / 2
			g_r_var = growth[i]*gamma_meas(r_var, w[1], w[2])/3600
		ELSEIF(r_back_ind[i+1]>r_nucl)
			// When r_var_1 surpasses value of r_nucl, set it equal to r_nucl
			// At this point, n(t_N) is calculated, used to obtain J.
			r_var_2 = r_back_ind[i+1]
			r_var_1 = r_nucl
			r_var = (r_var_1+r_var_2) / 2
			g_r_var = INTERP(r_nucl, r_back_ind, growth) /3600
			//j = i
		ELSE
			//To not let either r_var decrease below r_nucl.
			r_var_2 = r_nucl
			r_var_1 = r_nucl
			r_var = r_nucl
			g_r_var = interp(r_nucl, r_back_ind, growth) /3600
		ENDIF
		IF(g_r_var > 0)
			DF_wave[i] = (r_var_2/r_var_1)^(wall_const/g_r_var) // w[6] included in wall_const (see fit_Kw())
		ELSE
			DF_wave[i] = 1
		ENDIF
		
		// DF = deposition factor to correct n(a) backwards in time and size. 
		// DF = exp[-int(1/a)da*(wall_const/g_r_var)] = exp[ln(a1/a2)*(w_c/g)]
		// DF_wave[i] = 1/exp( -w[6]*wall_loss_diff(r_var)*(r_var_2-r_var_1)/g_r_var )
		// For coagulation, using average value of k(a) over delta_t (and thus delta_a)
		// integrate fitted power law for Kc(II) and sum over all values of A to obtain
		// pseudo first order rate constant, Kc1_r_var for Kc(a,A) for one value of a (r_var) and all values of A (radius).
		
		//Kc1_int1 = i_fit_c*r_var_1 + ( (c_fit_c/(p_fit_c+1)) * r_var_1^(p_fit_c+1) )
		//Kc1_int2 = i_fit_c*r_var_2 + ( (c_fit_c/(p_fit_c+1)) * r_var_2^(p_fit_c+1) )
		//Kc1_int = 0 //NaN
		//j = 0 // j = bin of n[j][i] that r_var is scavenged by.
		//DO // Only coag scavenging by particles larger than r_meas should be included.
		//	Kc1_int[j] = (n[j][i]+n[j][i+1])/2 * (Kc1_int2[j]-Kc1_int1[j]) * d_lnr[j]
		//	j = j+1
		//WHILE(radius[j] > r_meas)
		//Kc1_int = (n[p][i]+n[p][i+1])/2 * (Kc1_int2-Kc1_int1) * d_lnr
		//Kc2_r_var = w[3]*coag_rate_constant(r_var,radius) // 2nd order K_coag
		//Kc1 = (n[p][i]+n[p][i+1])/2 * Kc2_r_var * d_lnr	// ps. 1st order K_coag
		//Kc1_r_var = sum(Kc1_int, 0, L)
		//CF_wave[i] = exp( Kc1_r_var*(r_var_2-r_var_1)/g_r_var ) // for averaging method.
		//CF_wave[i] = exp( Kc1_r_var/g_r_var ) // for integration method.
		
		inner_int = 0					 // necessary to include before do- loop, to avoid Inf values.
		//inner_int = d_lnr * w[3]*coag_rate_constant(r_var,radius) * (n[p][i+1]+n[p][i])/2 // scav by all part.
		
		n_t1 = dNdlogDp[p][i]					// to read in from dN/dlnr matrix (produced by awk)
		n_t2 = dNdlogDp[p][i+1]				// dN/dlnr (units of r don't matter)
		replace_NaN_by_num(n_t1, 0)
		replace_NaN_by_num(n_t2, 0)
		//SMOOTH/B=2/F/E=3 7, n_t1	// NOTE THE SMOOTHING!!!   
		//SMOOTH/B=2/F/E=3 7, n_t2
	
		Smooth/B=2 5, n_t1
		Smooth/B=2 5, n_t2
		
		inner_int = d_lnr * w[3]*coag_rate_constant(r_var,radius) * (n_t1+n_t2)/2 		// scav by all part.
		
		//j = 0 							// j = bin of n[j][i] that r_var is scavenged by.
		//DO 							// include only coag scavenging by particles larger than r_meas, but not smaller than r_det
		//	inner_int[j] = d_lnr[j] * w[3]*coag_rate_constant(r_var,radius[j]) * (n[j][i+1]+n[j][i])/2
	 		//inner integral evaluated for r_var (average value of backwards calculated radius)
		//	j = j+1
		//WHILE( (j<L) && (radius[j]>r_var) ) 	//r_var = cohort radius; r_meas = measured radius
		////integrate inner_int_1 		// inner integral for r_var_1
		////integrate inner_int_2 		// inner integral for r_var_2
		INTEGRATE inner_int 		// inner integral for r_var
		IF(g_r_var > 0)
			CF_wave[i] = exp( inner_int[L-1] * (r_var_2 - r_var_1) / g_r_var )
		ELSE
			CF_wave[i] = 1
		ENDIF
		//print CF_wave[i]
		//print CF_wave[i-1]
		///////CF_wave[i] = exp( ( (inner_int_2[L-1]+inner_int_1[L-1])*(r_var_2-r_var_1) )/g_r_var )
		IF( r_back_ind[i+1]<r_nucl ) 	// if both radii are smaller than r_nucl,
			n_back_ind[i] = NaN 	// set n_back to NaN, to not disturb picture of n_back_t1
			//r_back_ind[i] = NaN
		ELSE
			n_back_ind[i] = n_back_ind[i+1] * ( DF_wave[i] * CF_wave[i] )
		ENDIF
		i = i-1
	WHILE((r_back_ind[i]>=r_nucl) && (i>=0) )
	//while( i>= 0) //54) //0 )
	// CF_wave and DF_wave are computed for each "backwards" timestep.
	// Total correction factor obtained by multiplying all individual ones.
	CF = 1
	DF = 1
	m = scan-1
	DO
		CF = CF * CF_wave[m]
		DF = DF * DF_wave[m]
		m = m-1
	WHILE(m>0)	
	//m = 0
	//DO
	//	CF = CF * CF_wave[m]
	//	DF = DF * DF_wave[m]
	//	m = m+1
	//WHILE(m<scan)	
	// loss is total correctioon factor for (pseudo) first order losses.
	// loss = exp[-(1/g)Int(k(a)da] for k = k_W + k_C(1) and integrated from r_nucl to r_meas.
	
	dep_F[bin][scan] = DF			// correction factor for removal by deposition (i.e. wall losses).
	coagsc_F[bin][scan] = CF		// correction factor for removal by coagulation scavenging.
	
	//wmcoag_F[bin][scan] = 1 + ( nn_meas*int_fb(r_meas, scan) ) // assumes steady state nucl.
	//wmcpulse_F[bin][scan] = 1+ ( N_total[scan]*int_f_pulse(r_meas, scan) ) // assumes nucl. pulse 
	// correction factor for within mode coagulation (removal and production)
	//b1_F[bin][scan] = calc_b1(r_nucl, r_meas, scan) // to compare with CF*DF
	//b1p_F[bin][scan] = calc_bp1(r_nucl, r_meas, scan) // to compare with CF*DF
	//return n_back_ind[number]/(CF*DF*radius[bin]) // dN/dR at t_form, first order corrected.
	//return n_back_ind[j]/radius[bin]		// dN/dR at , first order corrected.
	
	RETURN nn_meas*dep_F[bin][scan]*coagsc_F[bin][scan] 	//*wmcpulse_F[bin][scan]
	
	// dN/dR at t_form, corrected for coagulation and wall loss.
	// no interpolation necessary, since calculation is stopped once r_var exceeds 0.5
	//killwaves Kc1_int1, Kc1_int2, Kc1_int, Kc2_int, b, f_1, f_2, f
	//killwaves Kw_cor, Kc_cor, CF_wave, DF_wave, K_w_diff_fit, Kc2_r_var, Kc1, inner_int, inner_int_1, inner_int_2
END
  
 
  
FUNCTION fill_n_back(bin, scan) 	// Writes the output of the back calculation to matrices of radius and scan,
								// to enable reconstruction of size distributions as they are predicted to have existed.
	VARIABLE bin, scan
	NVAR K, L
	VARIABLE i
	WAVE n_back_ind, r_back_ind, radius, n_back_all, r_back_all
	n_back_all[(scan*L)+bin][] = n_back_ind[q]
	i = 0
	DO
		IF(r_back_ind[i] > 998)
			r_back_all[(scan*L)+bin][i] = NaN
		ELSE
			r_back_all[(scan*L)+bin][i] = r_back_ind[i]
		ENDIF
		i=i+1
	WHILE(i<K)
END


FUNCTION Get_Nucl_Rate_median(ctrlName) : ButtonControl      
 	STRING ctrlName
 	WAVE J_nucl, t_form_long
	Median_1Min_interval(J_nucl, t_form_long) 
END

FUNCTION Median_1Min_interval(J_nucl, t_form_long)	     
	WAVE t_form_long, J_nucl	
      	VARIABLE i, ref , m=0
      	WAVE t_wave
	MAKE/O/D/N=(5000) temp_wave
	
	print "Jnuc=", MedianEO(J_nucl, 0, DIMSIZE(t_form_long,0)-1)
	print "Jnuc=", Median(J_nucl, 0, DIMSIZE(t_form_long,0)-1)
	
	temp_wave[0]= t_wave[0]
	FOR(i=1;i<5000;i+=1)
		temp_wave[i] = temp_wave[i-1] + 60*10
	ENDFOR			 
	
	VARIABLE size = DIMSIZE(temp_wave, 0)   
      	VARIABLE size_inputdata = DIMSIZE(t_form_long, 0)   
        MAKE/O/D/N=(size) nucl_time_median, nucl_median
        nucl_time_median = NAN
        nucl_median = NAN  
       VARIABLE min_point, max_point   
        i=1 
        ref = 2    
    	DO 
        //IF( ref>1 )
        		min_point =  TRUNC( ( temp_wave[ref] + temp_wave[ref-1] ) / 2  ) 
           	max_point =  TRUNC( (temp_wave[ref] + temp_wave[ref+1] ) / 2  )    
          		
          	IF( t_form_long[i] > max_point  )          		           
          		DO
				ref = ref+1
          		        nucl_time_median [ref] = temp_wave[ref]
                                 nucl_median[ref] =   NAN 
          		        min_point =  TRUNC( ( temp_wave[ref] + temp_wave[ref-1] ) / 2  )           		 
          		 	max_point =  TRUNC( (temp_wave[ref] + temp_wave[ref+1] ) / 2  )
          		 WHILE(t_form_long[i] > max_point && ref < size)          		 
        		ENDIF

                IF(t_form_long[i] >= min_point   &&  t_form_long[i] <=max_point) 
                		MAKE/O/D temp_median
                        DO                                                                                      
				temp_median[m] = J_nucl [i]           
                                 m=m+1                                        
                                 i=i+1      
			WHILE( t_form_long[i] >= min_point   &&  t_form_long[i] <=max_point  && i< size_inputdata )
                         nucl_time_median [ref] = temp_wave[ref] 
                         IF(DIMSIZE(temp_median,0) >1)
                         	nucl_median[ref] =   MEAN( temp_median, 0, m-1) 
                         ELSE
                                  nucl_median[ref] = temp_median[0]       
                         ENDIF
                         KILLWAVES temp_median
                         m=0
		ELSEIF( t_form_long[i] > max_point  )
			ref = ref+1
		ELSEIF( t_form_long[i] < min_point)  
			i=i+1
		ENDIF                
	//ENDIF
        WHILE( i< size_inputdata  && ref < size)	 
        
        //Save/J/M="\r\n" nucl_median as "Jnucl_1min_median.txt"   
        //Save/J/M="\r\n" nucl_time_median as "Jnucl_1min_median_time.txt"   
END	
    




FUNCTION Get_Jmedian(ctrlName) : ButtonControl      
 	STRING ctrlName
 	NVAR K16
 	WAVE J_nucl, t_form_long
	Median_for_timeinterval(J_nucl, t_form_long) 
END

FUNCTION Median_for_timeinterval(J_nucl, t_form_long)	
			     
	WAVE t_form_long, J_nucl	
      	VARIABLE i, ref , m=0
      	WAVE t_wave  
      	NVAR K16
      	K16 = 1
	MAKE/O/D/N=(1500) temp_wave
	

	temp_wave[0]= t_wave[0]
	FOR(i=1;i<1500;i+=1)
		temp_wave[i] = temp_wave[i-1] + 60*K16
	ENDFOR			 
	
	VARIABLE size = DIMSIZE(temp_wave, 0)   
      	VARIABLE size_inputdata = DIMSIZE(t_form_long, 0)   
        MAKE/O/D/N=(size) new_nucl_time_median, new_nucl_median
        new_nucl_time_median = NAN
        new_nucl_median = NAN 
       VARIABLE min_point, max_point
        i=1 
        ref = 2    
    	DO 
        //IF( ref>1 )
        		min_point =  TRUNC( ( temp_wave[ref] + temp_wave[ref-1] ) / 2  ) 
           	max_point =  TRUNC( (temp_wave[ref] + temp_wave[ref+1] ) / 2  )    
          		
          	IF( t_form_long[i] > max_point  )          		           
          		DO
				ref = ref+1
          		        new_nucl_time_median [ref] = temp_wave[ref]
                                 new_nucl_median[ref] =   NAN 
          		        min_point =  TRUNC( ( temp_wave[ref] + temp_wave[ref-1] ) / 2  )           		 
          		 	max_point =  TRUNC( (temp_wave[ref] + temp_wave[ref+1] ) / 2  )
          		 WHILE(t_form_long[i] > max_point && ref < size)          		 
        		ENDIF

                IF(t_form_long[i] >= min_point   &&  t_form_long[i] <=max_point) 
                		MAKE/O/D temp_median
                        DO                                                                                      
				temp_median[m] = J_nucl [i]           
                                 m=m+1                                        
                                 i=i+1      
			WHILE( t_form_long[i] >= min_point   &&  t_form_long[i] <=max_point  && i< size_inputdata )
                         new_nucl_time_median [ref] = temp_wave[ref] 
                         IF(DIMSIZE(temp_median,0) >1)
                         	new_nucl_median[ref] =   MEAN( temp_median, 0, m-1) 
                         ELSE
                                  new_nucl_median[ref] = temp_median[0]       
                         ENDIF
                         KILLWAVES temp_median
                         m=0
		ELSEIF( t_form_long[i] > max_point  )
			ref = ref+1
		ELSEIF( t_form_long[i] < min_point)  
			i=i+1 
		ENDIF                
	//ENDIF
        WHILE( i< size_inputdata  && ref < size)	    
END	



FUNCTION Nucl_Rate_Graphs(ctrlName) : ButtonControl      
 	STRING ctrlName 	
 	variable x1 
	J_nucl_graph()	
END

FUNCTION J_nucl_graph() : Graph
	PauseUpdate; Silent 1		// building window...
	
	WAVE growth,t_midgrowth
	DUPLICATE/D growth growth_new
	VARIABLE i
	FOR(i=0;i<DIMSIZE(t_midgrowth,0);i=i+1)
		IF(growth[i] < -300 || growth[i] > 300 )  
			growth_new[i] = NAN
		ELSE
			growth_new[i] = growth[i]
		ENDIF  
	ENDFOR     
	DUPLICATE/D growth_new growth_smth
	SMOOTH/B=3 3, growth_smth   
	
	Display growth_new vs t_midgrowth
	SetAxis bottom 36000,75600
	SetAxis left 0,40
	ModifyGraph mode=3,marker(growth_new)=5,rgb(growth_new)=(0,0,0)
	AppendToGraph/R J_nucl vs t_form_long
	ModifyGraph mode(J_nucl)=3
	ModifyGraph msize(J_nucl)=2.5,rgb(J_nucl)=(30464,30464,30464)
	AppendToGraph growth_smth vs t_midgrowth
	ModifyGraph mode(growth_smth)=4,marker(growth_smth)=16,lsize(growth_smth)=1;DelayUpdate
	ModifyGraph rgb(growth_smth)=(0,0,65280)
	
	AppendToGraph nucl_median vs nucl_time_median
	ModifyGraph mode(nucl_median)=4,marker(nucl_median)=7;DelayUpdate
	ModifyGraph rgb(nucl_median)=(65280,0,0)
	SetAxis right 0.001,100  
	Legend/C/N=text0/F=0/A=MC
	
	ModifyGraph hideTrace(growth_new)=1,hideTrace(J_nucl)=1
	Label left "GR (nm h-1)"
	Label right "Jnuc (cm-3 s-1)"   
	ModifyGraph hideTrace=0
	
END


 






END