#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

//Plot Bivariate contour ellipse for arbitrary xy data
//John Economides 06-06-19

Function plotBCE(xData,yData,CI) //CI is intended Confidence Interval from 0 to 1

Wave/Z xData,ydata
Variable CI
Variable ellipseArea
Variable majorAxis,minorAxis
Variable psiDeg

Display yData vs xData

WaveStats/Q xData
Variable xAvg=v_avg
Variable xSD=v_sdev
WaveStats/Q yData
Variable yAvg=v_avg
Variable ySD=v_sdev

Variable aspectatio=(Wavemax(yData)-Wavemin(yData))/(Wavemax(xData)-Wavemin(xData))

Variable pearsons=StatsCorrelation(xData,yData)
Variable covariance=pearsons*xSD*ySD

Make/O/N=(2,2)/D covMat
covMat[0][0]=xSD^2
covMat[1][1]=ySD^2
covMat[0][1]=covariance
covMat[1][0]=covariance

MatrixEigenV covMat
Wave W_eigenValues
Variable a=sqrt(real(W_eigenValues[0]))
Variable b=sqrt(real(W_eigenValues[1]))

Variable psi=.5*(atan((1/aspectatio)*((2*covariance)/((xSD^2)-(ySD^2)))))
Variable srInvChis=sqrt(StatsInvChiCDF(CI,2))

Make/O/N=50/D theta,ellipseX,ellipseY
theta=2*pi*p/49
ellipseY=((b*cos(psi)*sin(theta))+(a*cos(theta)*sin(psi)))*srInvChis+yAvg
ellipseX=((a*cos(theta)*cos(psi))-(b*sin(theta)*sin(psi)))*srInvChis+xAvg

AppendToGraph ellipseY vs ellipseX
ModifyGraph mode=3, marker=19, msize=1
ModifyGraph rgb=(0,0,0,25000)

ModifyGraph mode(ellipseY)=0
ModifyGraph rgb(ellipseY)=(65000,0,0)

ModifyGraph width={Plan,1,bottom,left}

ellipseArea=pi*(a*srInvChis)*(b*srInvChis)
majorAxis=a*srInvChis
minorAxis=b*srInvChis
psiDeg=psi*180/pi

Print "major axis length = " + num2str(majorAxis)
Print "minor axis length = " + num2str(minorAxis)

if (Wavemax(ellipseY)-Wavemin(ellipseY)>Wavemax(ellipseX)-Wavemin(ellipseX))
	Print "angle psi = " + num2str(psiDeg+90)
else
	Print "angle psi = " + num2str(psiDeg)
endif

Print "ellipse area = " + num2Str(ellipseArea)


End