#pragma rtGlobals=1		// Use modern global access method.
#include "massdist"
#include "HIPPO_cloudcut"
#include "insertnans"
#include "meanbin"
//run cloud cutting program first to generate cloud_flag wave.
function one_second_final_data(flightdate, analysisdate)
	STRING flightdate, analysisdate // in format "MM DD YYYY"
	// this function was written to do a complete final-data analysis of SP2 data for ARCPAC. 
	//Modified for preliminary data submission for Macpex.
	// it requires:
	// the official AOCTIMEWAVE and ambient temperature and static pressure.  These should be read from the 
	// most current aircraft met file on the archive:
	// aoctimewave, s since midnight of start date
	// ambient_temp: deg C   
	// static_press: mBarr
	// also...
	// SP2 data in the format below, on a single timebase called "time_s" from the .mhk file:
	// sample_dp_1: vccs
	// sample_t_1:  deg C
	// sample_p_1:  torr
	// Yag_power_1, in volts. This guy is used only to remove data with the laser turned off: do not smooth
	// also..
	// SP2 data from raw data processed in labview:
	// bstart_time, s: on the aoctimewave, but missing entries. perhaps short, perhaps long. 
	// binc_mass:  in femtograms, binned to AOCTIMEWAVE, on bstart_time  
	// bduty_cycle; as a fraction of 1, binned to aoctimewave.  This is just the number of 0.2 s buffer reads that contributed mass to each binc_mass 
	//                                                                                           divided by 5. 
	//
	// This function will do the rest:
	
	
	//Rewritten 2012/01/10 for HIPPO IV+V analysis using aircract data stream from GV
	wave bstart_time, time_s, ambient_temp, static_press, time_hms, aoctime
	
	//get rid of NaNs in ambient temp and static_press
	duplicate /o ambient_temp ambt
	duplicate /o time_s ts
	ts = numtype(ambient_temp[p])==0 ? time_s[p]:NaN
	wavetransform zapnans ts
	wavetransform zapnans ambt
	
	ambient_temp = interp(time_s,ts, ambt)
	
	duplicate /o static_press ambt
	duplicate /o time_s ts
	ts = numtype(static_Press[p])==0 ? time_s[p]:NaN
	wavetransform zapnans ts
	wavetransform zapnans ambt
	
	static_press = interp(time_s,ts, ambt)

	// Declare waves:
	wave sample_p_1, sample_dp_1, binc_mass, sample_t_1, yag_power_1
	wave aoctime, time_s, bduty_cycle,bstart_time 
	make /o/n=(numpnts(sample_p_1)) smth_times = NaN
	
	//Find cloud flag. If it doesn't exist, make one that won't cut anything but warn the user.
	if (exists("cloud_flag")==1)
		wave cloud_flag
	else 
	print "Cloud_flag wave not found! Did you run the cloud_cut analysis? Making a dummy vector of zeros."
		make /n=(numpnts(time_s)) cloud_flag=0
	endif
	
	duplicate /o bstart_time sample_dp_terp
	sample_dp_terp = interp(bstart_time,time_s, sample_dp_1)
	
	variable /g n=0
	variable /g j=0
	
	do

		if (sample_dp_terp[n]<=0)
			binc_mass[n] = NaN
			binc_mass[n+1] = NaN
			j+=1
		endif
		n=n+1
	
	while (n<numpnts(bstart_time))
	
	print "removed "+num2str(j)+" points for negative flow"
	variable jstor = j
	
	// Calculate timeshift to AOCTIMEWAVE, and clean up time_s:
	//find_sp2_aoc_timeshift()  // the timeshift is called: sp2_aoc_time_shift
	find_sp2_aoc_timeshift_hms()
	variable /g sp2_aoc_time_shift
	duplicate /o aoctime, aoctimewave
	print 1
	
	// First, get all the binned BC data/dutycycle information on the AOCTIMEWAVE scale:
	make/o/n=(numpnts(aoctimewave)) inc_mass_temp = nan, duty_cycle = nan
	variable i = 0

	do  // get rid of changes of day in the start_time:
		if(bstart_time[i+1] < bstart_time[i])
			bstart_time[i+1] += 24*3600
		endif
		i+=1
	while(i<numpnts(bstart_time)-1)
	print 2
	i = 0
	do                     // insert missing points into the LV processed sp2 data. 
		if(bstart_time[i+1] > bstart_time[i]+1) // missed the next point:
			insertpoints i+1,1, bstart_time, binc_mass, bduty_cycle
			bstart_time[i+1] = bstart_time[i] +1  // now i is the new, inserted point
			binc_mass[i+1] = nan
			bduty_cycle[i+1] = nan
			j+=1
		endif
		
		i+=1
	while(i<numpnts(bstart_time)-1)
		
	print num2str(j-jstor)+" NaNs inserted for missing bstart_times"
	jstor = j
		
	variable /g n = 0 ,m = 0          // after this bit of code, aoctimewave[n] and bstart_time[m] will represent the earliest time at which we can start filling
	if(aoctimewave[0] < bstart_time[0]) // in real data.
		do
			n+=1
			j+=1
		while(aoctimewave[n] < bstart_time[0]) 
	else
		do
			m+=1
		while(aoctimewave[0] > bstart_time[m])
	endif 

print 3
print num2str(j-jstor)+" NaNs inserted for bstart_time[0] > aoctimewave[0]"
jstor = j

	do  // after this loop all the usable data has been put into the aoctimewave space.  If SP2 wasn't working, the code above put in
		inc_mass_temp[n] = binc_mass[m]   // nans to take care of it.  If Sp2 kept running after landing, it doesn't matter. 
		duty_cycle[n] = bduty_cycle[m]

		n+=1 
		m+=1
	while(n< numpnts(aoctimewave) && m < numpnts(bstart_time))

print 4
	// Now, smooth the HK data streams, and convert to correct units:
	wave smth
	smooth_wave(sample_dp_1, time_s, 3)
	duplicate/o smth sample_dp_smth
	smooth_wave(sample_p_1, time_s, 3)
	duplicate/o smth sample_p_smth
	smooth_wave(sample_t_1, time_s,3)
	duplicate/o smth sample_t_smth
	sample_t_smth += 273.15  // convert to K
	Yag_power_1 = abs(yag_power_1)
print 5
	// the smoothed data time wave is called "smth_times", Shift it so that it is in aoctimewave space:
	smth_times -= sp2_aoc_time_shift  + 0.0 // for 415 flight: add 0.9 sec offset to co2
	/////////////////////////////////////////////
	// CRITICAL!!!!!    CRITICAL!!!!!!!!  CRITICAL!!!!!!!!!!!/
	/////////////////////////////////////////////
	// above: adjust smth_times to reflect the time delay in transport of aerosol to the SP2.  I have looked at 
	// our data of 4/18 flight against CO2, and determined that there was no additional time delay beyond that 
	// determined from the aoc/sp2 time shift.  That suggests that we're lucky and the time delay of AOC data on the
	// network matched the aerosol delay at our position behind the LTI!  However, in general, this is not the case:  
	// for 4.12 flight: add 1 for total of 3.5
	// for 415 flight: add 0.9 sec offset to co2 for total offset of 3.5
	// for 4.18 flight add 0 for total offset of 3.15
	//for 4.19 flight add 0 for total offset of 3.9
	// for 4.21 flight add0 for total of 2.28
	// 4.23 flight: add 1 for total of 4.6

	// Now, interpolate these values to AOCTIMEWAVE.  Note that if SP2 wasn't running, the interpolation gives crap, but as inc_mass is nan, it doesn't matter.
	make/o/n=(numpnts(aoctimewave)) sample_flow_vccs, sample_p_mbarr, sample_t_k, ambient_t_k, staticprs
	make/o/n=(numpnts(aoctimewave))  yag_power_v, cloud_flag_terp
	sample_flow_vccs = interp(aoctimewave, smth_times, sample_dp_smth)
	sample_p_mbarr = interp(aoctimewave, smth_times, sample_p_smth)*1013/760 //Converted sample p into Torr.
	sample_t_k = interp(aoctimewave, smth_times, sample_t_smth)
	ambient_t_k =  interp(aoctimewave, time_s, ambient_temp)
	ambient_t_k+=273.15
	staticprs =  interp(aoctimewave, time_s, static_press)
	yag_power_v = interp(aoctimewave, time_s , yag_power_1)
	cloud_flag_terp = interp(aoctimewave, time_s, cloud_flag)
	
	if (exists("masscorr")==2)
		variable /g masscorr
	else
		print "Couldn't find mass correction factor. Applying no correction."
		variable /g masscorr = 1
	endif

	// generate our final data product:
	make/o/n=(numpnts(aoctimewave)) bc_ng_kg= nan, bc_ng_m3 = nan

	// note that the binc_mass is already corrected for the n of m writing issues 
	// timestep of 1 sec is left out:
	bc_ng_m3 = inc_mass_temp*sample_t_k*staticprs/sample_flow_vccs/ambient_t_k/sample_p_mbarr/duty_cycle*masscorr
	bc_ng_kg = inc_mass_temp*sample_t_k*1013/sample_flow_vccs/273.15/sample_p_mbarr/1.29/duty_cycle*masscorr

	// Null values where:
	//1) aoctimewave out of SP2 operating range. Put in above
	//2) sample_flow_vccs is below 0.4 vccs
	//3) clouds are present
	//4) Yag_power_v is less than 1
	// 5) SP2 was not running  within flight: this is nulled by LV programs.
	// 6) PMT1 less than some amount: I did the translation to aoc wave manually like this:
	//////duplicate aoctimewave pmt1_onaoc
	//////„duplicate time_s time_s_shifted
	//////„time_s_shifted = time_s - sp2_aoc_time_shift
	//////„pmt1_onaoc = interp(aoctimewave, time_s_shifted, pmt_read1)
//	wave pmt1_onaoc
	n =0

	do 

		//2:
		if(sample_flow_vccs[n] < 0.4 || sample_flow_vccs[n] > 5.5)
			bc_ng_kg[n] = nan
			bc_ng_m3[n] = nan
			j+=1
		endif
		
		//2a.
		if (sample_dp_terp[n]<=0)
			bc_ng_kg[n] = NaN
			bc_ng_kg[n+1] = NaN
			
			bc_ng_m3[n] = NaN
			bc_ng_m3[n+1] = NaN
			j+=2
		endif

		// 3: 
		//if(droplet_flag[n] > 0.1)
			//bc_ng_kg[n] = nan
			//bc_ng_m3[n] = nan
		//endif


		//4: 
		if(yag_power_v[n] <0.2)
			bc_ng_kg[n] = nan
			bc_ng_m3[n] = nan
			j+=1
		endif

		// 5. 
//		if(pmt1_onaoc[n] < 5.285)
	//		bc_ng_kg[n] = nan
		//	bc_ng_m3[n] = nan
//		endif

		//6: 
		if(cloud_flag_terp[n] > 0.5)
			bc_ng_kg[n] = nan
			bc_ng_m3[n] = nan
			j+=1
		endif
		n+=1


	while(n < numpnts(aoctimewave))

	
	print "Removed "+num2str(j-jstor)+" points for sample flow, yag power and clouds"
	jstor = j
	print "Total points removed = "+num2str(jstor)
	variable /g shifttime
	aoctimewave -= shifttime
	// now print the archive file:
	HIPPOSP2file(flightdate, analysisdate)

end





Function HIPPOSP2file(Flightdate,Analysisdate)
	String Flightdate,Analysisdate
	String /g flightmonth = Flightdate[0,1], analysismonth = Analysisdate[0,1] 
	String /g flightday = Flightdate[3,4], analysisday = Analysisdate[3,4]
	String /g flightyear = Flightdate[6,9], analysisyear = Analysisdate[6,9]
	wave aoctimewave, bc_ng_kg, bc_ng_m3
	//Prompt fdatestring, "Enter flight and process date string (yyyy mm dd yyyy mm dd): " 
	Variable spfile, n, m							//make refNum variable
	open /M="Enter archive filename: " spfile		//create and open file (0verwrite ifit exists)
	fprintf spfile, "26	1001\n"     					// # of lines in header (NLHEAD) and file type (1001)
	fprintf spfile, "Perring, A.E. (anne.perring@noaa.gov),  Schwarz, J. P., Spackman, J.R., Gao, R.-S., Fahey, D. W., Watts, L. A.\n"
	fprintf spfile, "NOAA Earth System Research Laboratory, Chemical Sciences Division, Boulder, Colorado\n"
	fprintf spfile, "NOAA Single-Particle Soot Photometer (SP2)\n"     // instrument name (SNAME)
	fprintf spfile, "HIPPO-5, NCAR/NSF G-V, 09 August - 08 September 2011\n"     					// mission name (MNAME)
	fprintf spfile, "1	1 \n"								//File volume number and number of volumes
	fprintf spfile, "%s	%s	%s	%s	%s	%s \n" FLightyear, Flightmonth, Flightday, ANalysisyear, Analysismonth, Analysisday
	fprintf spfile, "0 \n"
	fprintf spfile, "UTC, seconds since midnight on flight day \n"
	fprintf spfile, "2 \n"     							// number of dependent variables
	fprintf spfile, "1	1 \n"     						// scale factor
	fprintf spfile, "-9999.99	-9999.99\n"    			// NAN variable
	fprintf spfile, "BC_ng_kg, Black carbon mass mixing ratio: ng-BC/kg-air \n"
	fprintf spfile, "BC_ng_m3, Black carbon volume mixing ratio: ng-BC/m3 of air at ambient conditions \n"
	fprintf spfile, "1 \n"  
	fprintf spfile, "10 \n"  
	fprintf spfile, "REVISION: R1\n"  
	fprintf spfile, "R1: Final 1 second data. \n"  
	fprintf spfile, "*************************************************   Final Data   **********************************************\n"   	// special comment lines
	fprintf spfile, "The SP2 detects refractory black carbon (BC) in the accumulation mode (~90 - 600 nm) covering about 90 pct of \n"    							 
	fprintf spfile, "the mass distribution. The data here have been scaled by a factor of 1.08 to represent the total BC loading in the \n"
	fprintf spfile, "accumulation mode. BC data acquired in clouds have been removed based on SP2 internal diagnostics, the 2D-C and \n"
	fprintf spfile, " CDP cloud probe, the liquid water sensor (PLWC), the raw icing-rate indicator (RICE), relative humidity and \n"
	fprintf spfile, "dew point temperature measurements. The measurement uncertainty associated with the data is 30 pct. \n"
	fprintf spfile, "**************************************************************************************************************** \n"
	fprintf spfile, "UTC	BC_ng_kg	BC_ng_m3 \n"
		
	// Convert NaN to -9999
	n=0
	do 
		if (numtype(bc_ng_kg[n]) != 0  %^ numtype(bc_ng_m3[n]) !=0)
			bc_ng_kg[n] = -9999.99
			bc_ng_m3[n]= -9999.99
		endif

		n+=1
	while (n< numpnts(aoctimewave))
	
	wfprintf spfile, "%6.0f	%7.2f	%7.2f \n" aoctimewave, bc_ng_kg, bc_ng_m3
	//two space delimited integers and two floats with automatic width and 3 digits precision
	close spfile
	
End



function find_sp2_aoc_timeshift()
	wave time_s, time_hms
	// first, generate the aoctimewave to compare to sp2 time (i.e. time_s):
	duplicate/o time_hms, aoctime

	// correct for roll-over days
	variable n = 0, lastaoc
	do
		if (numtype(aoctime[n]) == 0)
			// real number
			lastaoc = aoctime[n]
		endif
		if(aoctime[n+1] < lastaoc)
			aoctime[n+1] +=86400
		endif
		n+=1
	while(n < numpnts(aoctime)-1)

	// check for rollover in time_s:
	n = 0
	do
		if(time_s[n+1] - time_s[n]<0)
			time_s[n+1] += 24*3600
		endif 
		n+=1 
	while(n<numpnts(time_s)-1)

	// find the average timeshift -- if time_s is quite stably measured every second, 
	duplicate/o aoctime tmp               //then this shift could be off  by upto 1 second. 
	tmp= time_s - aoctime
	wavestats/q tmp
	print "The average timeshift (SP2 - AOC) is : ", v_avg, " seconds"
	variable/g sp2_aoc_time_shift = v_avg
end

//Added 1/10/12 because new housekeeping file reads "Time_HMS" variable as seconds rather than as hour minute and second.
function find_sp2_aoc_timeshift_hms()
wave time_hms, time_s
wavetransform zapnans time_hms
duplicate /o time_hms, aoctime
// correct for roll-over days
	variable n = 0, lastaoc
	do
		if (numtype(aoctime[n]) == 0)
			// real number
			lastaoc = aoctime[n]
		endif
		if(aoctime[n+1] < lastaoc)
			aoctime[n+1] +=86400
		endif
		n+=1
	while(n < numpnts(aoctime)-1)

	// check for rollover in time_s:
	n = 0
	do
		if(time_s[n+1] - time_s[n]<0)
			time_s[n+1] += 24*3600
		endif 
		n+=1 
	while(n<numpnts(time_s)-1)

	// find the average timeshift -- if time_s is quite stably measured every second, 
	duplicate/o aoctime tmp               //then this shift could be off  by upto 1 second. 
	tmp= time_s - aoctime
	wavestats/q tmp
	print "The average timeshift (SP2 - AOC) is : ", v_avg, " seconds"
	variable/g sp2_aoc_time_shift = v_avg

end

Function Smooth_wave(wa,times,step)
	// you give wave wa on time base times, specify steps, and get smth and smth_times
	// step is the number to average together
	wave wa, times
	variable step 
	variable n = 0 ,m=0,cave,tave,num=0, hol

	make/o/n= (floor(numpnts(wa)/step)) smth =nan, smth_times =nan
	do
		cave = 0
		tave = 0
		m=0
		hol = 0
		do
	
			hol +=1
			cave += wa[n]
			tave += times[n]
	
			m+=1
			n+=1		
		while(m<step)

		smth[num] = cave/hol
		smth_times[num] = tave/hol
		num+=1
	while(n< numpnts(wa))
end

i=numpnts(timewave_new)
do 
destwave[i]=NaN
i-=1
while (timewave_new[i]>timewave_orig[inf])

end