#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

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/G a=sqrt(real(W_eigenValues[0]))
Variable/G 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)

Print "ellipse area = " + num2Str(ellipseArea)

End