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


Menu "Macros"
	"Get ScaleFactor Between Waves",GetScaleFactorBetweenWaves(0,0,"","","","")
	"re-bin wave",rebin_wave("",10,1)
End


Function rebin_wave(inwave,nchan,renormalize)
	String inwave					// name of wave to bin up
	Variable nchan					// number of points to bin together
	Variable renormalize			// set to 0 or 1
	if (!(nchan>0) || numtype(renormalize) || exists(inwave)!=1)
		nchan = (!(nchan>0)) ? 10 : nchan
		renormalize = numtype(renormalize) ? 0 : renormalize
		inwave = SelectString(exists(inwave)==1,"",inwave)
		Prompt inwave,"input wave for rebinning",popup, WaveList("*", ";", "")
		Prompt nchan,"number of points in input wave to sum together"
		Prompt renormalize,"renormalize result",popup,"true;false"
		renormalize += 1
		DoPrompt "re-bin wave",inwave,nchan,renormalize
		if (V_flag || nchan<=1)
			return 1
		endif
		renormalize -= 1
	endif

	Wave in = $inwave
	String outwave = inwave+"_"+num2istr(nchan)
	Variable nin=numpnts(in)
	Variable nout = ceil(nin/nchan)			//number of points in outwave wave
	Variable lo_edge						//used for scaling the $outwave
	String W_xunits = StringByKey("XUNITS",WaveInfo($inwave, 0))

	Duplicate/O $inwave $outwave
	SetScale/P x 0,1,"", $outwave
	Wave out = $outwave
	out = sum(out,p*nchan,p*nchan+nchan-1)
	Redimension/N=(nout) $outwave
	lo_edge = (nchan-1)/2 * deltax(in) + leftx(in)	//find left edge for scaling
	SetScale/P x, (lo_edge), (deltax(in)*nchan), W_xunits, $outwave
	if (renormalize)
		out /= nchan						//renormalize the binned wave
	endif
End



//Macro rebin_wave(inwave,nchan,renormalize)
//	String inwave
//	Prompt inwave,"input wave for rebinning",popup, WaveList("*", ";", "")
//	Variable nchan=10
//	Prompt nchan,"number of points in input wave to sum together"
//	Variable renormalize=2
//	Prompt renormalize,"renormalize result",popup,"true;false"
//
//	Silent 1 ;  DelayUpdate
//	if (nchan<=1)
//		return
//	endif
//
//	String outwave = inwave+"_"+num2istr(nchan)
//	Variable xx,ii
//	Variable nin=numpnts($inwave)
//	Variable nout = ceil(nin/nchan)			//number of points in outwave wave
//	Variable lo_edge						//used for scaling the $outwave
//
//	Variable i1,i2							//string indicies
//	String Wunits = WaveInfo($inwave, 0)	//obtain both units for the SetScale
//	i1 = strsearch(Wunits,"XUNITS:",0) + 7
//	i2 = strsearch(Wunits, ";",i1) - 1
//	String W_xunits = Wunits[i1,i2]
//
//	duplicate/O $inwave $outwave
//	SetScale/P x 0,1,"", $outwave
//	$outwave = sum($outwave,p*nchan,p*nchan+nchan-1)
//	Redimension/N=(nout) $outwave
//	lo_edge = (nchan-1)/2 * deltax($inwave) + leftx($inwave)	//find left edge for scaling
//	SetScale/P x, (lo_edge), (deltax($inwave)*nchan), W_xunits, $outwave
//	if (renormalize==1)
//		$outwave /= nchan								//renormalize $outwave
//	endif
//End Macro


Function GetScaleFactorBetweenWaves(i1,i2,yy1,xx1,yy2,xx2)
	Variable i1, i2			// range in points of (yy1,xx1) to use
	String xx1,yy1			// first (x,y) pair
	String xx2,yy2			// second (x,y) pair
	Prompt i1, "first point to use of Wave1"
	Prompt i2, "last point to use of Wave1"
	Prompt yy1,"Wave1 Y values",popup,WaveList("*",";","")
	Prompt xx1,"Wave1 X values",popup,WaveList("*",";","")
	Prompt yy2,"Wave2 Y values",popup,WaveList("*",";","")
	Prompt xx2,"Wave2 X values",popup,WaveList("*",";","")
	DoPrompt "get scale factor",i1,i2,yy1,xx1,yy2,xx2
	if (V_flag)
		return 1
	endif
	if (!monotonic($xx2))
		Abort "the wave "+xx2+" is NOT monotonic"
	endif
	Variable scale = ScaleFactorBetweenWaves(i1,i2,$yy1,$xx1,$yy2,$xx2)
	print "  the scaling is:       "+yy1+" = "+yy2+" *",scale
End
//Macro GetScaleFactorBetweenWaves(i1,i2,cc1,xx1,cc2,xx2)
//	Variable i1, i2			// range in points of (cc1,xx1) to use
//	Prompt i1, "first point to use of Wave1"
//	Prompt i2, "last point to use of Wave1"
//	String xx1,cc1			// first (x,y) pair
//	Prompt cc1,"Wave1 Y values"
//	Prompt xx1,"Wave1 X values"
//	String xx2,cc2			// second (x,y) pair
//	Prompt cc2,"Wave2 Y values"
//	Prompt xx2,"Wave2 X values"
//	Silent 1
//	if (!monotonic($xx2))
//		Abort "the wave "+xx2+" is NOT monotonic"
//	endif
//	Variable scale = ScaleFactorBetweenWaves(i1,i2,$cc1,$xx1,$cc2,$xx2)
//	print "  the scaling is:       "+cc1+" = "+cc2+" *",scale
//EndMacro
Function ScaleFactorBetweenWaves(i1,i2,cc1,xx1,cc2,xx2)
	Variable i1, i2			// range in points of (cc1,xx1) to use
	Wave xx1,cc1			// first (x,y) pair
	Wave xx2,cc2			// second (x,y) pair
	//
	//	This routine based on min of   SUM{ [Yi - a*f(Xi)]^2 } = V
	//	dV/da = SUM{ 2*[Yi - a*f(Xi)] * [-f(Xi) }
	//	dV/da = -2 * SUM{ Yi*f(Xi) - a*f(Xi)^2 }
	//	for dV/da = 0 implies
	//	0 = -2 * SUM{ Yi*f(Xi) - a*f(Xi)^2 }
	//	0 = SUM{ Yi*f(Xi) - a*f(Xi)^2 }
	//	0 = SUM{ Yi*f(Xi) } - a*SUM{ f(Xi)^2 }
	//	SUM{ Yi*f(Xi) } = a*SUM{ f(Xi)^2 }
	//
	//	a = SUM{ Yi*f(Xi) } / SUM{ f(Xi)^2 }
	//
	Variable yf		// SUM{ Yi*f(Xi) }
	Variable f2		// SUM{ f(Xi)^2 }
	Variable Ncc1 = numpnts(cc1)
	if (!WaveExists(cc1) || !WaveExists(xx1) || !WaveExists(cc2) || !WaveExists(xx2))
		return NaN
	endif

	if ((i2<1) || (i2>=Ncc1))
		i2 = Ncc1-1
	endif

	Variable fXi		// f(Xi)
	yf = 0
	f2 = 0
	Variable i=i1
	do
		fXi = interp(xx1[i], xx2, cc2)
		yf += cc1[i] * fXi
		f2 += fXi^2
		i += 1
	while (i<=i2)
	return yf/f2
End




// added proper handling of NaN's in ywave
Function CombineNearbyPoints(yName,xName,yErr,dx,firstX)
	String yName			// ="Det"
	String xName			// ="dE"
	String yErr				// ="Det_err"
	Variable dx				// max difference between two points considered the same (=0.05)
	Variable firstX			// =-1e30

	Prompt yName,"Name of Y wave",popup, WaveList("*", ";", "")
	Prompt xName,"Name of X wave",popup, WaveList("*", ";", "")
	Prompt yErr,"Name of Y error wave",popup, "_none_;_root_N_;"+WaveList("*", ";", "")
	Prompt dx, "Max difference between two points considered the same"
	Prompt firstX, "X value of where to start Averaging"

	dx = abs(dx)
	if (exists(yName)!=1)
		yName=""
	endif
	if (exists(xName)!=1)
		xName=""
	endif
	if (strlen(yName)<1 || strlen(xName)<1 || dx<=0 || numtype(firstX)==2)
		DoPrompt "Combine nearby points", yName,xName,yErr,dx,firstX
	endif
	PauseUpdate

	Variable err = (cmpstr("_none_",yErr)) && (strlen(yErr)>0)		// true when yErr wave given
	Variable root_N=(!cmpstr("_root_N_",yErr))						// used sqrt(N) for error
	if (root_N)
		err = 0
	endif

	if (exists(yName)!=1 || exists(xName)!=1)
		Abort "trying to combine waves that do not exist!"
	endif
	Wave yw=$yName
	Wave xw=$xName
	if (numpnts(yw)!=numpnts(xw))
		Abort "'"+yName+"' and '"+xName+"'  must be the same length"
	endif
	if (err)
		Wave yew=$yErr
		if (numpnts(yw)!=numpnts(yew))
			Abort "'"+yName+"' and '"+yErr+"'  must be the same length"
		endif
	endif
	if (!cmpstr(yName,xName))
		Abort "xwave and ywave must be different"
	endif

	if (err)
		sort xw,xw,yw,yew
	else
		sort xw,xw,yw
	endif

	// remove all NaN's in either xw, yw, or yew
	Variable i = numpnts(xw)-1								// start checking from the end of the waves
	if (err)
		do
			if (numtype(xw[i]+yw[i]+yew[i])==2)// found a NaN
				DeletePoints i, 1, yw,xw,yew
			endif
			i -= 1
		while(i>=0)
	else
		do
			if (numtype(xw[i]+yw[i])==2)		// found a NaN
				DeletePoints i, 1, yw,xw
			endif
			i -= 1
		while(i>=0)
	endif

	Variable sumY,sumX,n,j
	Variable sumE, sigma2
	i=max(0,BinarySearch(xw,firstX))
	do
		if( abs(xw[i+1]-xw[i])<dx )
			n = 0
			do					// make n the number of points to combine it should be at least 2
				n += 1
			while( (abs(xw[i+n]-xw[i])<dx) && ((i+n)<numpnts(xw)))
			sumY = 0
			sumX = 0
			sumE = 0
			sigma2 = 1
			j = 0
			do					// sum up n points and reduce length of xName, yName, and yErr
				if (err)
					sigma2 = (yew[j+i])^2
				endif
				if (root_N)
					sigma2 = abs(yw[j+i])			// sigma = sqrt(y[i])
				endif
				sumX += xw[j+i]
				sumY += yw[j+i] / sigma2
				sumE += 1/sigma2
				j += 1
			while (j<n)

			xw[i] = sumX/n
			xw[(i+1),] = xw[p+n-1]
			yw[i] = sumY / sumE
			yw[(i+1),] = yw[p+n-1]
			if (err)
				yew[i] = sqrt(1/sumE)
				yew[(i+1),] = yew[p+n-1]
				Redimension/N=(numpnts(xw)-n+1) yw,xw,yew
			else
				Redimension/N=(numpnts(xw)-n+1) yw,xw
			endif
		endif

		i += 1
	while (i<(numpnts(xw)-1))
End
Function CombineNearbyPointsOLD(yName,xName,yErr,dx,firstX)
	String yName			// ="Det"
	String xName			// ="dE"
	String yErr				// ="Det_err"
	Variable dx				// max difference between two points considered the same (=0.05)
	Variable firstX			// =-1e30

	Prompt yName,"Name of Y wave",popup, WaveList("*", ";", "")
	Prompt xName,"Name of X wave",popup, WaveList("*", ";", "")
	Prompt yErr,"Name of Y error wave",popup, "_none_;_root_N_;"+WaveList("*", ";", "")
	Prompt dx, "Max difference between two points considered the same"
	Prompt firstX, "X value of where to start Averaging"

	dx = abs(dx)
	if (exists(yName)!=1)
		yName=""
	endif
	if (exists(xName)!=1)
		xName=""
	endif
	if (strlen(yName)<1 || strlen(xName)<1 || dx<=0 || numtype(firstX))
		DoPrompt "Combine nearby points", yName,xName,yErr,dx,firstX
	endif
	PauseUpdate

	Variable err = (cmpstr("_none_",yErr)) && (strlen(yErr)>0)		// true when yErr wave given
	Variable root_N=(!cmpstr("_root_N_",yErr))						// used sqrt(N) for error
	if (root_N)
		err = 0
	endif

	if (exists(yName)!=1 || exists(xName)!=1)
		Abort "trying to combine waves that do not exist!"
	endif
	Wave yw=$yName
	Wave xw=$xName
	if (numpnts(yw)!=numpnts(xw))
		Abort "'"+yName+"' and '"+xName+"'  must be the same length"
	endif
	if (err)
		Wave yew=$yErr
		if (numpnts(yw)!=numpnts(yew))
			Abort "'"+yName+"' and '"+yErr+"'  must be the same length"
		endif
	endif
	if (!cmpstr(yName,xName))
		Abort "xwave and ywave must be different"
	endif

	if (err)
		sort xw,xw,yw,yew
	else
		sort xw,xw,yw
	endif

	Variable Ndups=0			// number of duplicate sets found
	Variable sumY,sumX,n,j,i=BinarySearch(xw,firstX)
	Variable sumE, sigma2
	i = max(0,i)
//	print "starting at point",i
	do

		if( abs(xw[i+1]-xw[i])<dx )
			Ndups += 1
			n = 0
			do					// make n the number of points to combine it should be at least 2
				n += 1
			while( (abs(xw[i+n]-xw[i])<dx) && ((i+n)<numpnts(xw)))
			sumY = 0
			sumX = 0
			sumE = 0
			sigma2 = 1
			j = 0
			do					// sum up n points and reduce length of xName, yName, and yErr
				if (err)
					sigma2 = (yew[j+i])^2
				endif
				if (root_N)
					sigma2 = abs(yw[j+i])			// sigma = sqrt(y[i])
				endif
				sumX += xw[j+i]
				sumY += yw[j+i] / sigma2
				sumE += 1/sigma2
				j += 1
			while (j<n)

			xw[i] = sumX/n
			xw[(i+1),] = xw[p+n-1]
			yw[i] = sumY / sumE
			yw[(i+1),] = yw[p+n-1]
			if (err)
				yew[i] = sqrt(1/sumE)
				yew[(i+1),] = yew[p+n-1]
				Redimension/N=(numpnts(xw)-n+1) yw,xw,yew
			else
				Redimension/N=(numpnts(xw)-n+1) yw,xw
			endif
		endif

		i += 1
	while (i<(numpnts(xw)-1))
//	print "    found ",Ndups," duplicate sets"
End




//  use				 #include <Concatenate Waves>


//// Concatenate Waves version 1.01
//// 		Feb. 2, 1997, JW added test for text waves
////		Nov. 20, 2000, JZT modified to use StringFromList & NumberByKey, rather than <Strings as Lists>
//
////	ConcatenateWaves(w1, w2)
////		Tacks the contents of w2 on the to end of w1.
////		If w1 does not exist, it is created.
////		This is designed for 1D waves only.
//Function ConcatenateWaves(w1, w2)
//	String w1, w2
//	
//	Variable numPoints1, numPoints2
//
//	if (Exists(w1) == 0)
//		Duplicate $w2, $w1
//	else
//		String wInfo=WaveInfo($w2, 0)
//		Variable WType=NumberByKey("NUMTYPE", wInfo)
//		numPoints1 = numpnts($w1)
//		numPoints2 = numpnts($w2)
//		Redimension/N=(numPoints1 + numPoints2) $w1
//		if (WType)				// Numeric wave
//			Wave/C ww1=$w1
//			Wave/C ww2=$w2
//			ww1[numPoints1, ] = ww2[p-numPoints1]
//		else						// Text wave
//			Wave/T tw1=$w1
//			Wave/T tw2=$w2
//			tw1[numPoints1, ] = tw2[p-numPoints1]
//		endif
//	endif
//End
//
////	ConcatenateWavesInList(dest, wl)
////		Makes a dest wave that is the concatenation of the source waves.
////		Overwrites the dest wave if it already exists.
////		wl is assumed to contain at least one wave name.
////		This is designed for 1D waves only.
//Function ConcatenateWavesInList(dest, wl)
//	String dest		// name of output wave
//	String wl		// semicolon separated list of waves ("w0;w1;w2;")
//	
//	Variable i					// for walking through wavelist
//	String theWaveName
//	Variable destExisted
//	
//	destExisted = Exists(dest)
//	if (destExisted)
//		Redimension/N=0 $dest
//	endif
//	
//	i = 0
//	do
//		theWaveName = StringFromList(i, wl, ";")
//		if (strlen(theWaveName) == 0)
//			break										// no more waves
//		endif
//		if (cmpstr(theWaveName, dest) != 0)		// don't concat dest wave with itself
//			ConcatenateWaves(dest, theWaveName)
//		endif
//		i += 1
//	while (1)
//End



// this revision of gcf uses the the "PrimeFactors" command available with Igor 5
// is also assumes that 0 can be divided by anything, so gcf of {0,4,8} is 4
Function gcf(factors)				// find greatest common factor of a wave whose values are all integers
	String factors					// either a list of factors or the name of a wave containing the factors

	String tempName=""
	Variable N,i
	if (exists(factors)==1 && strsearch(factors,";",0)<0)		// passed a wave name
		Wave wav=$factors
	else
		N = ItemsInList(factors)								// a list was passed
		if (N<1)
			return NaN
		endif
		tempName = UniqueName("tempPrimes",1,0)
		Make/N=(N) $tempName
		Wave wav=$tempName
		for (i=0;i<N;i+=1)
			wav[i] = str2num(StringFromList(i,factors))
		endfor
	endif
	N=numpnts(wav)
	if (N<1)
		return NaN
	endif

	String primeListsName = UniqueName("ListOfPrimes",1,0)
	Make/N=(N)/T $primeListsName
	Wave/T primeLists=$primeListsName
	primeLists = ""
	String allPrimes="", item
	Variable j,gcf
	for (i=0;i<N;i+=1)
		PrimeFactors/Q  wav[i]
		Wave W_PrimeFactors=W_PrimeFactors
		for (j=0;j<numpnts(W_PrimeFactors);j+=1)
			primeLists[i]=AddListItem(num2istr(W_PrimeFactors[j]), primeLists[i])
		endfor
		allPrimes+= primeLists[i]
		if (wav[i]==0)									// ensure "0" is there for a 0
			primeLists[i]=AddListItem("0", primeLists[i])
		endif
	endfor

	Variable m,found,jN = ItemsInList(allPrimes)
	j = 0
	gcf=1
	for (j=0;j<jN;j+=1)								// loop over all of the primes
		item = StringFromList(j,allPrimes)
		found = 1
		for (i=0;i<N;i+=1)
			m = WhichListItem(item,primeLists[i])
			if (m>=0)
				primeLists[i] = RemoveListItem(m,primeLists[i])
			elseif (cmpstr(primeLists[i],"0;"))
				found=0
			endif
		endfor
		if (found)
			gcf *= str2num(item)
		endif
	endfor
	KillWaves/Z primeLists,W_PrimeFactors, $tempName
	return gcf
End
//// this revision of gcf uses the the "PrimeFactors" command available with Igor 5
//Function gcf(factors)				// find greatest common factor of a wave whose values are all integers
//	String factors					// either a list of factors or the name of a wave containing the factors
//
//	String tempName=""
//	Variable N,i
//	if (exists(factors)==1 && strsearch(factors,";",0)<0)		// passed a wave name
//		Wave wav=$factors
//	else
//		N = ItemsInList(factors)
//		if (N<1)
//			return NaN
//		endif
//		tempName = UniqueName("tempPrimes",1,0)
//		Make/N=(N) $tempName
//		Wave wav=$tempName
//		for (i=0;i<N;i+=1)
//			wav[i] = str2num(StringFromList(i,factors))
//		endfor
//	endif
//
//	Variable j
//	N=numpnts(wav)
//	if (N<1)
//		return NaN
//	endif
//	String primeListsName = UniqueName("ListOfPrimes",1,0)
//	Make/N=(N)/T $primeListsName
//	Wave/T primeLists=$primeListsName
//	primeLists = ""
//	String allPrimes="", item
//	for (i=0;i<N;i+=1)
//		PrimeFactors/Q  wav[i]
//		Wave W_PrimeFactors=W_PrimeFactors
//		for (j=0;j<numpnts(W_PrimeFactors);j+=1)
//			primeLists[i]=AddListItem(num2istr(W_PrimeFactors[j]), primeLists[i])
//		endfor
//		allPrimes+= primeLists[i]
//	endfor
//
//	Variable m,found,jN = ItemsInList(allPrimes),gcf=1
//	for (j=0;j<jN;j+=1)
//		item = StringFromList(j,allPrimes)
//		found = 1
//		for (i=0;i<N;i+=1)
//			m = WhichListItem(item,primeLists[i])
//			if (m>=0)
//				primeLists[i] = RemoveListItem(m,primeLists[i])
//			else
//				found=0
//			endif
//		endfor
//		if (found)
//			gcf *= str2num(item)
//		endif
//	endfor
//	KillWaves/Z primeLists,W_PrimeFactors, $tempName
//	return gcf
//End
//
//
//// this routine only works on integres that fit in unsigned 32bit ints
//Function gcf(wav)				// find greatest common factor of a wave whose values are all integers
//	Wave wav
//
//	Variable N=numpnts(wav)
//	if (N<1)
//		return NaN
//	endif
//
//	String intsName1 = UniqueName("ints",1,0)
//	Make/N=(N)/U/I $intsName1
//	String intsName2 = UniqueName("ints",1,0)
//	Make/N=(N)/U/I $intsName2
//	Wave ints = $intsName1
//	Wave remain = $intsName2
//	ints = round(abs(wav[p]))					// convert from input form to integres
//	WaveStats/Q ints
//	if (V_max<1)
//		KillWaves/Z $intsName1,$intsName2
//		return 1
//	endif
//
//	MakeListOfPrimes(V_max)
//	Wave ListOfPrimes=root:ListOfPrimes
//	Variable nPrimes=numpnts(ListOfPrimes)
//	Variable tolerance = 0.1/V_max
//	Variable factor = 1
//	Variable prime, ip
//	for (ip=0;ip<nPrimes;ip+=1)
//		prime = ListOfPrimes[ip]				// get next prime
//		do										// remove as many factors of prime as possible
//			remain = mod(ints[p],prime)
//			if (sum(remain,0,N-1)<tolerance)	// all of the values are divisibl by prime
//				ints /= prime
//				factor *= prime
//			else
//				break
//			endif
//		while(1)
//	endfor
//	KillWaves/Z $intsName1,$intsName2
//	return factor
//End
// test(0,0,0)
// test(21,-49,7)
// test(21,21,7)
// test(20,21,7)
//Function test(a,b,c)
//	Variable a,b,c
//	Make/N=3/O test3_test
//	test3_test={a,b,c}
//	print gcf(test3_test), gcf3(a,b,c)
//End
//
//Static Function MakeListOfPrimes(N)
//	Variable N		// N can be up to 65536, but no greater
//	if (N>65536)
//		Abort "can only make list of primes up to 65536"
//	endif
//	if (exists("root:ListOfPrimes")==1)		// if list exists, check its length
//		if (N<=NumberByKey("max",note($"root:ListOfPrimes"),"="))
//			return numpnts($"root:ListOfPrimes")		// existing list is long enough, do not re-compute
//		endif
//	endif
//
//	Make/N=6542/O/U/I root:ListOfPrimes
//	Wave list = root:ListOfPrimes
//	list = NaN
//	Note/K list
//	Note list, "max="+num2istr(N)
//
//	Variable i,j=0
//	for (i=2;i<=N;i+=1)
//		if (TestForPrime(i))
//			list[j] = i
//			j += 1
//		endif
//	endfor
//	Redimension/N=(j) list
//	return j
//End
//Static Function TestForPrime(n)
//	Variable n
//	Variable nmax = ceil(sqrt(n))
//	n = round(abs(n))
//	if (n<2.5)
//		return 1
//	endif
//	Variable i
//	for (i=2;i<=nmax;i+=1)
//		if (mod(n,i)<1e-8)		// found a divisor
//			return 0
//		endif
//	endfor
//	return 1
//End