#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3				// Use modern global access method and strict wave access
#pragma DefaultTab={3,20,4}		// Set default tab width in Igor Pro 9 and later

// Procedures to generate 3D plots of real H wave functions orbitals
//	  Melissa A. Hines, Cornell University, Melissa.Hines@cornell.edu

Menu "Macros"
	"Generate Bohr Orbital"
	"Generate H Orbital"
	"Modify Isosurface"
End

Proc GenerateHOrbital(n, l, m, genPointCloud)
variable n, l, m, genPointCloud
Prompt n, "n:"
Prompt l, "l:"
Prompt m, "m:"
Prompt genPointCloud, "Generate Point Cloud?", popup "Yes;No"
	
	makeWavefunctions(0, 0, n, l, m, genPointCloud==1)
End

Proc GenerateBohrOrbital(n, amplitude)
variable n, amplitude=0.2

	make/n=256/o circle_x, circle_y, bohrOrbital_x, bohrOrbital_y; 
	SetScale/I x 0,2*pi,"", circle_x,circle_y,bohrOrbital_x,bohrOrbital_y
	circle_x = cos(x);circle_y = sin(x)
	bohrOrbital_x = (1+amplitude*cos(n*x))*cos(x);bohrOrbital_y =(1+amplitude*cos(n*x))*sin(x)
	
	if(strlen(WinList("bohrOrbitalGraph", ";","")) < 2)
    	Execute "bohrOrbitalGraph()"
    endif


End

Proc ModifyIsosurface(theSurface, larger, factor)
variable theSurface, larger, factor=2
Prompt theSurface, "Which isosurface?", popup "Psi;Psi2"
Prompt larger, "Which direction?", popup "Larger;Smaller"
Prompt factor, "How much?"

	PauseUpdate
	variable oldIsoValue, multiplier
	
	multiplier = (larger < 2) ? 1/factor : factor
	
	if(theSurface < 2)	// Modify Psi
		oldIsoValue = currentIsoValue("GizmoPsi", "posSurface")
	
		DoWindow/F GizmoPsi
    	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={isoValue,oldIsoValue*multiplier}
    	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={isoValue,oldIsoValue*multiplier}
    else
		oldIsoValue = currentIsoValue("GizmoPsi2", "psi2Surface")
	
		DoWindow/F GizmoPsi2
    	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={isoValue,oldIsoValue*multiplier}
    endif
End

// Calculate H atom wavefunction using real orbitals defined at
//    https://en.wikipedia.org/wiki/Atomic_orbital#Real_orbitals
Function calcPsi(n, l, m, rr, theta, phi)
variable n, l, m, rr, theta, phi

    if(m < 0)
    	return sqrt(2)*(-1)^m * Rnl(n, l, rr) * imag(sphericalHarmonics(l,-m,theta,phi))
    elseif(m > 0)
    	return sqrt(2)*(-1)^m * Rnl(n, l, rr) * real(sphericalHarmonics(l,m,theta,phi))
    else
    	return Rnl(n, l, rr) * real(sphericalHarmonics(l,0,theta,phi))
    endif
End

// Calculate H atom radial wavefunction as defined at https://en.wikipedia.org/wiki/Hydrogen_atom#:~:text=The%20normalized%20position%20wavefunctions
Function Rnl(n, l, r)
variable n, l, r

	variable out, a0 = 0.529, rho
	
	rho = 2*r/n/a0
	out = sqrt((2/n/a0)^3*factorial(n - l - 1)/(2 * n * factorial(n + l))) * exp(-rho/2)*(rho)^(l) * LaguerreA(n-l-1, 2*l+1, rho)
	return out 
End

// This function calculates the radial wavefunction and displays graphs of psi, psi^2, and 4 pi r^2 psi^2.
// Three dimensional waves representing psi^2 and positive and negative regions of psi are generated. These
//   are displayed in Gizmo plots. The isosurfaces may need tweaking to generate a pleasing size.
// Finally (if requested) a Monte Carlo algorithm is used to generate a point cloud representing psi.
Function makeWavefunctions(pts, range, n, l, m, doMakeCloud)
// pts = number of points in sqr box. 75 seems to be a good number. Use 0 to autoscale.
// range = Box is defined over the region x, y, z = (-range, range). Use 0 to autoscale.
// n, l, m: Usual hydrogenic quantum numbers
// doMakeCloud = (0, 1) depending on whether the point cloud should be generated.
variable pts, range, n, l, m, doMakeCloud
 
    variable i,j,k,xx,yy,zz,rr,rho,phi,theta, psi
    
    // Validate the quantum numbers
    if( (n < 1) || (l >= n) || (l < 0) || (abs(m) > l))
    	Print "Invalid quantum numbers"
    	return -1
    endif
    
    // Autosizing if pts == 0
    if(pts < 1)
    	pts = 75
    endif
    variable center=floor(pts/2)	// Use odd number of points
    
    // Autoscaling if range == 0
    if(!exists("radialPsi"))
    	make/n=512/o radialPsi, radialPsi2, radialProbability, zero    
    endif
    wave radialPsi, radialPsi2, radialProbability, zero
    if(range < 0.01)    	
    	range = calcRange(150, n, l)	// Set a huge range and do a rough graph
    	range = calcRange(range, n, l)	// Refine the range
    endif  
    
    // Update radial wavefunction
    SetScale/I x 0,range,"", radialPsi,radialPsi2,radialProbability,zero
    radialPsi = Rnl(n,l,x)
    radialPsi2 = radialPsi^2
    radialProbability = 4*pi*x^2*radialPsi2
    
    // If the radial graphs are not displayed, display them
    if(strlen(WinList("radialPsiGraph", ";","")) < 2)
    	Execute "radialPsiGraph()"
    endif
    if(strlen(WinList("radialPsi2Graph", ";","")) < 2)
    	Execute "radialPsi2Graph()"
    endif
    if(strlen(WinList("radialProbabilityGraph", ";","")) < 2)
    	Execute "radialProbabilityGraph()"
    endif
    	

 	pts = 2*floor(pts/2) + 1
    
    make/o/n=(pts, pts, pts) ddd_psi2, ddd_psiPos, ddd_psiNeg
    SetScale/I x -range,range,"", ddd_psi2, ddd_psiPos, ddd_psiNeg
	SetScale/I y -range,range,"", ddd_psi2, ddd_psiPos, ddd_psiNeg
	SetScale/I z -range,range,"", ddd_psi2, ddd_psiPos, ddd_psiNeg 
    for(k=0;k<pts;k+=1)
        zz=range * (k-center)/center 
        for(j=0;j<pts;j+=1)
            yy=range * (j-center)/center
            for(i=0;i<pts;i+=1)
                xx=range * (i-center)/center
                
                // compute spherical polar coords
                rr=sqrt(xx^2 + yy^2 + zz^2)
                phi=atan2(yy,xx)
                theta=acos(zz/rr)
                
                psi = calcPsi(n, l, m, rr, theta, phi)
                if(psi > 0)
                	ddd_psiPos[i][j][k] = psi
                	ddd_psiNeg[i][j][k] = 0
                else
                	ddd_psiPos[i][j][k] = 0
                	ddd_psiNeg[i][j][k] = -psi
                endif
                ddd_psi2[i][j][k] = psi^2
            endfor
        endfor
    endfor
    
    if(strlen(WinList("GizmoPsi", ";","")) < 2)
    	Execute "GizmoPsi()"
    endif
    if(strlen(WinList("GizmoPsi2", ";","")) < 2)
    	Execute "GizmoPsi2()"
    endif

    wavestats/q ddd_psiPos
    DoWindow/F GizmoPsi
    ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={isoValue,V_max/5}
    ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={isoValue,V_max/5}

    wavestats/q ddd_psi2
    DoWindow/F GizmoPsi2
    ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={isoValue,V_max/50}
    
    if(doMakeCloud)
	    if(strlen(WinList("GizmoPointCloud", ";","")) > 2)
    		DoWindow/F GizmoPointCloud
	    	ModifyGizmo stopUpdates
   		endif	    
	    makeCloud(n, l, m, range, 30000)
	    if(strlen(WinList("GizmoPointCloud", ";","")) > 2)
    		ModifyGizmo resumeUpdates
    	else 
    		Execute "GizmoPointCloud()"
   		endif
	endif
	
	wave W_coef, W_sigma
	KillWaves/Z W_coef, W_sigma, scaledWave, M_colors 
End

// Calculates an optimal range for the wavefunctions based on the outermost maximum in the radial distribution function
Function calcRange(roughRange, n, l)
variable roughRange, n, l

    wave radialPsi, radialPsi2, radialProbability, zero

	// Set a huge range and do a rough graph
	SetScale/I x 0,roughRange,"", radialPsi,radialPsi2,radialProbability,zero
    radialPsi = Rnl(n,l,x)
    radialPsi2 = radialPsi^2
    radialProbability = 4*pi*x^2*radialPsi2
    
    // Fit the max in the radialProbabily to a Gaussian
    wavestats/Q radialProbability
    K0 = 0;K1 = V_max;K2 = V_maxLoc;
    CurveFit/H="1110"/Q gauss radialProbability
    wave W_coef
    return ceil(W_coef[2] + 1.5 * W_coef[3])
End

// Returns the currently plotted isoValue for a Gizmo isosurface
Function currentIsoValue(gizmoName, isoSurfaceName)
string gizmoName, isoSurfaceName

	string matchStr, theStr
	variable cur = 0
	
	GetGizmo/n=$gizmoName objectList
	wave/T TW_gizmoObjectList
	
	matchStr = "*ModifyGizmo ModifyObject="+isoSurfaceName+",objectType=isoSurface,property={ isoValue,*"	
	do
		if(stringmatch(TW_gizmoObjectList[cur], matchStr))
			theStr = TW_gizmoObjectList[cur]
			theStr = theStr[strsearch(theStr, "isoValue,", 0)+9,inf]
			theStr = theStr[0, strlen(theStr)-2]
			KillWaves TW_gizmoObjectList
			return str2num(theStr)
		endif
		cur += 1
	while(cur < numpnts(TW_gizmoObjectList))
	KillWaves TW_gizmoObjectList
	return NaN	
End

Function colorXYZ(n, l, m, xx, yy, zz)
variable n, l, m, xx, yy, zz

	variable rr, phi, theta
	rr=sqrt(xx^2 + yy^2 + zz^2)
    phi=atan2(yy,xx)
    theta=acos(zz/rr)
	return sign(calcPsi(n, l, m, rr, theta, phi))

End

Function acceptOnePnt(n, l, m, xx, yy, zz, numAccepted)
variable n, l, m, xx, yy, zz, numAccepted

	wave pointCloud, colorPnts
	pointCloud[numAccepted][0] = xx; pointCloud[numAccepted][1] = yy; pointCloud[numAccepted][2] = zz;
	colorPnts[numAccepted] = colorXYZ(n, l, m, xx, yy, zz)
	numAccepted += 1
	return numAccepted
End

Function acceptPnts(n, l, m, xx, yy, zz, numAccepted)
variable n, l, m, xx, yy, zz, numAccepted

	// Take advatage of the symmetry of the orbitals to calculate faster
	numAccepted = acceptOnePnt(n, l, m, xx, yy, zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, -xx, yy, zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, xx, -yy, zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, xx, yy, -zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, -xx, -yy, zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, -xx, yy, -zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, -xx, -yy, -zz, numAccepted)
	numAccepted = acceptOnePnt(n, l, m, -xx, -yy, -zz, numAccepted)
	return numAccepted
End

// This function uses a Monte Carlo algorithm to generate a probabalistic point cloud representing the wavefunction.
// The positive and negative components of the wave function are colored red and blue, respectively.
Function makeCloud(n, l, m, range, numAccepts)
variable n, l, m, range, numAccepts
// range = Box is initiallydefined over the region x, y, z = (-range, range). Autoscales after 10% complete

	variable xx, yy, zz, rr, phi, theta, numTries, numAccepted,magnitude, kMax, maxRange, rescale
	numAccepts = 8 * ceil(numAccepts/8)
	make/n=(numAccepts, 3)/o pointCloud
	pointCloud = 0
	make/o/n=(numAccepts) colorPnts
	numAccepted = 0
	numTries = 0
	rescale = 1
	
	// Find the maximum value of the wavefuction over the selected range
	kMax = findMax(n, l, m, numAccepts, range)
	
	// The efficiency of the algorithm depends on an optimal setting of the range variable.
	// We use the input range for 10% of the tries, then reset based on the max value
	// found to that point
	maxRange = -1
	do
		xx = enoise(range); yy = enoise(range); zz = enoise(range)
		rr=sqrt(xx^2 + yy^2 + zz^2)
        phi=atan2(yy,xx)
        theta=acos(zz/rr)
		magnitude = calcPsi(n, l, m, rr, theta, phi)^2
		if(abs(magnitude) > abs(enoise(kMax)))
			numAccepted = acceptPnts(n, l, m, xx, yy, zz, numAccepted)
			if(rescale && max(xx, yy, zz) > maxRange)
				maxRange = max(xx, yy, zz)
			endif
		endif
		numTries += 1
		if(rescale && numAccepted > 0.1 * numAccepts)
			range = maxRange * 1.15
			rescale = 0
			print "New range =", range
		endif
	while(numAccepted < numAccepts)
	
	// To create a color wave for a scatter in Gizmo
	ColorTab2Wave Rainbow
	wave M_colors
	duplicate/o colorPnts, scaledWave
	Variable theMin=WaveMin(scaledWave)
	Variable theMax=WaveMax(scaledWave)
	Variable nor=(DimSize(M_colors,0)-1)/2
	MatrixOP/O scaledWave=nor*(scaledWave-theMin)
	M_colors/=65535
	
	make/o/n=(numAccepts,4) colorWave=1 // alpha will be 1.
	colorWave[][0]= M_colors[scaledWave[p]][0]
	colorWave[][1]= M_colors[scaledWave[p]][1]
	colorWave[][2]= M_colors[scaledWave[p]][2]
	
	// If calculating radialPsi^2, the following line should be approx 1
	//print "The integral is ", numAccepted/numTries * 8 * range^3 * kMax
	// print numTries, numAccepted

End

Function findMax(n, l, m, numTries, range)
variable n, l, m, numTries, range

	variable theMax = 0, cnt = 0, xx, yy, zz, rr, phi, theta, magnitude
	do
		xx = enoise(range); yy = enoise(range); zz = enoise(range)
		rr=sqrt(xx^2 + yy^2 + zz^2)
        phi=atan2(yy,xx)
        theta=acos(zz/rr)
        
		magnitude = abs(calcPsi(n, l, m, rr, theta, phi))
		if(magnitude > theMax)
			theMax = magnitude
		endif
		cnt += 1
	while(cnt < numTries)

	return theMax
End
Window GizmoPsi() : GizmoPlot
	PauseUpdate; Silent 1		// building window...
	// Building Gizmo 9 window...
	NewGizmo/W=(35,40,550,500)
	ModifyGizmo startRecMacro=901
	ModifyGizmo scalingOption=63
	AppendToGizmo isoSurface=root:ddd_psiPos,name=posSurface
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ surfaceColorType,0}
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ lineColorType,0}
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ lineWidthType,0}
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ fillMode,2}
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ lineWidth,1}
	ModifyGizmo ModifyObject=posSurface,objectType=isoSurface,property={ isoValue,0.01}
	ModifyGizmo modifyObject=posSurface,objectType=Surface,property={calcNormals,1}
		ModifyGizmo setObjectAttribute={posSurface,diffusePos}
		ModifyGizmo setObjectAttribute={posSurface,shininess0}
		ModifyGizmo setObjectAttribute={posSurface,specular}
	AppendToGizmo Axes=boxAxes,name=axes0
	ModifyGizmo ModifyObject=axes0,objectType=Axes,property={-1,axisScalingMode,1}
	ModifyGizmo ModifyObject=axes0,objectType=Axes,property={0,ticks,3}
	ModifyGizmo ModifyObject=axes0,objectType=Axes,property={1,ticks,3}
	ModifyGizmo ModifyObject=axes0,objectType=Axes,property={2,ticks,3}
	ModifyGizmo modifyObject=axes0,objectType=Axes,property={-1,Clipped,0}
	AppendToGizmo isoSurface=root:ddd_psiNeg,name=negSurface
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ surfaceColorType,0}
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ lineColorType,0}
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ lineWidthType,0}
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ fillMode,2}
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ lineWidth,1}
	ModifyGizmo ModifyObject=negSurface,objectType=isoSurface,property={ isoValue,0.01}
	ModifyGizmo modifyObject=negSurface,objectType=Surface,property={calcNormals,1}
		ModifyGizmo setObjectAttribute={negSurface,diffuseNeg}
		ModifyGizmo setObjectAttribute={negSurface,shininess0}
		ModifyGizmo setObjectAttribute={negSurface,specular}
	AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
	ModifyGizmo modifyObject=freeAxesCue0,objectType=freeAxisCue,property={colorType,0}
	AppendToGizmo light=Directional,name=light0
	ModifyGizmo modifyObject=light0,objectType=light,property={ position,-0.184655,-0.414741,0.891007,0.000000}
	ModifyGizmo modifyObject=light0,objectType=light,property={ direction,-0.184655,-0.414741,0.891007}
	ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
	AppendToGizmo attribute diffuse={1,0,0,1,1032},name=diffusePos
	AppendToGizmo attribute diffuse={1.5259e-05,0.244434,1,1,1032},name=diffuseNeg
	AppendToGizmo attribute shininess={42,42},name=shininess0
	AppendToGizmo attribute specular={1,1,1,1,1032},name=specular
	ModifyGizmo setDisplayList=0, object=freeAxesCue0
	ModifyGizmo setDisplayList=1, object=light0
	ModifyGizmo setDisplayList=2, object=posSurface
	ModifyGizmo setDisplayList=3, object=negSurface
	ModifyGizmo autoscaling=1
	ModifyGizmo currentGroupObject=""
	ShowTools
	ModifyGizmo showInfo
	ModifyGizmo infoWindow={551,31,1185,418}
	ModifyGizmo SETQUATERNION={0.000000,-0.000000,0.000000,1.000000}
	ModifyGizmo endRecMacro
	Execute/Q/Z "SetWindow kwTopWin sizeLimit={46,234,inf,inf}" // sizeLimit requires Igor 7 or later
EndMacro

Window GizmoPsi2() : GizmoPlot
	PauseUpdate; Silent 1		// building window...
	// Building Gizmo 9 window...
	NewGizmo/W=(213,605,728,1065)
	ModifyGizmo startRecMacro=901
	ModifyGizmo scalingOption=63
	AppendToGizmo isoSurface=root:ddd_psi2,name=psi2Surface
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ surfaceColorType,0}
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ lineColorType,0}
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ lineWidthType,0}
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ fillMode,2}
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ lineWidth,1}
	ModifyGizmo ModifyObject=psi2Surface,objectType=isoSurface,property={ isoValue,0.000225374}
	ModifyGizmo modifyObject=psi2Surface,objectType=Surface,property={calcNormals,1}
		ModifyGizmo setObjectAttribute={psi2Surface,diffuse}
		ModifyGizmo setObjectAttribute={psi2Surface,specular}
		ModifyGizmo setObjectAttribute={psi2Surface,shininess}
	AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
	ModifyGizmo modifyObject=freeAxesCue0,objectType=freeAxisCue,property={colorType,0}
	AppendToGizmo light=Directional,name=light0
	ModifyGizmo modifyObject=light0,objectType=light,property={ position,0.294889,0.694714,0.656059,0.000000}
	ModifyGizmo modifyObject=light0,objectType=light,property={ direction,0.294889,0.694714,0.656059}
	ModifyGizmo modifyObject=light0,objectType=light,property={ specular,1.000000,1.000000,1.000000,1.000000}
	AppendToGizmo attribute diffuse={1,0.250019,0.250019,1,1032},name=diffuse
	AppendToGizmo attribute specular={1,1,0,1,1032},name=specular
	AppendToGizmo attribute shininess={30,30},name=shininess
	ModifyGizmo setDisplayList=0, object=freeAxesCue0
	ModifyGizmo setDisplayList=1, object=light0
	ModifyGizmo setDisplayList=2, object=psi2Surface
	ModifyGizmo setDisplayList=3, opName=clearColor0, operation=clearColor, data={1,1,1,1}
	ModifyGizmo autoscaling=1
	ModifyGizmo currentGroupObject=""
	ShowTools
	ModifyGizmo showInfo
	ModifyGizmo infoWindow={782,750,1405,1094}
	ModifyGizmo SETQUATERNION={0.000000,0.000000,0.000000,1.000000}
	ModifyGizmo endRecMacro
	Execute/Q/Z "SetWindow kwTopWin sizeLimit={46,234,inf,inf}" // sizeLimit requires Igor 7 or later
EndMacro

Window GizmoPointCloud() : GizmoPlot
	PauseUpdate; Silent 1		// building window...
	// Building Gizmo 9 window...
	NewGizmo/W=(947,721,1512,1219)
	ModifyGizmo startRecMacro=901
	ModifyGizmo scalingOption=63
	AppendToGizmo Scatter=root:pointCloud,name=thePoints
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ scatterColorType,1}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ markerType,0}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ sizeType,0}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ rotationType,0}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ Shape,2}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ size,0.05}
	ModifyGizmo ModifyObject=thePoints,objectType=scatter,property={ colorWave,root:colorWave}
	AppendToGizmo freeAxesCue={0,0,0,1},name=freeAxesCue0
	ModifyGizmo modifyObject=freeAxesCue0,objectType=freeAxisCue,property={colorType,0}
	ModifyGizmo setDisplayList=0, object=freeAxesCue0
	ModifyGizmo setDisplayList=1, object=thePoints
	ModifyGizmo autoscaling=1
	ModifyGizmo currentGroupObject=""
	ShowTools
	ModifyGizmo showInfo
	ModifyGizmo infoWindow={804,539,1398,783}
	ModifyGizmo SETQUATERNION={-0.080336,0.084591,0.010790,0.993114}
	ModifyGizmo endRecMacro
	Execute/Q/Z "SetWindow kwTopWin sizeLimit={46,234,inf,inf}" // sizeLimit requires Igor 7 or later
EndMacro

Window radialPsiGraph() : Graph
	PauseUpdate; Silent 1		// building window...
	Display /W=(4,523,399,731) radialPsi,zero
	ModifyGraph gFont="Helvetica"
	ModifyGraph lStyle(zero)=1
	ModifyGraph rgb=(0,0,0)
	ModifyGraph mirror=2
	ModifyGraph nticks(left)=0
	ModifyGraph font="Helvetica"
	ModifyGraph fSize=12
	ModifyGraph lblMargin(left)=11
	ModifyGraph lblLatPos(left)=-1
	ModifyGraph btLen=4
	Label left "ψ"
	Label bottom "\\f02r\\f00 (Å)"
EndMacro

Window radialPsi2Graph() : Graph
	PauseUpdate; Silent 1		// building window...
	Display /W=(0,771,395,979) radialPsi2
	ModifyGraph gFont="Helvetica"
	ModifyGraph rgb=(0,0,0)
	ModifyGraph mirror=2
	ModifyGraph nticks(left)=0
	ModifyGraph font="Helvetica"
	ModifyGraph fSize=12
	ModifyGraph lblMargin(left)=10
	ModifyGraph btLen=4
	Label left "ψ\\S2\\M"
	Label bottom "\\f02r\\f00 (Å)"
EndMacro

Window radialProbabilityGraph() : Graph
	PauseUpdate; Silent 1		// building window...
	Display /W=(1,1015,396,1223) radialProbability
	ModifyGraph gFont="Helvetica"
	ModifyGraph rgb=(0,0,0)
	ModifyGraph mirror=2
	ModifyGraph nticks(left)=0
	ModifyGraph font="Helvetica"
	ModifyGraph fSize=14
	ModifyGraph lblMargin(left)=11
	ModifyGraph btLen=4
	Label left "4 π \\f02r\\f00\\S2\\M ψ\\S2\\M"
	Label bottom "\\f02r\\f00 (Å)"
EndMacro

Window bohrOrbitalGraph() : Graph
	PauseUpdate; Silent 1		// building window...
	Display /W=(35,520,424,903) circle_y vs circle_x
	AppendToGraph bohrOrbital_y vs bohrOrbital_x
	ModifyGraph gFont="Helvetica",height={Plan,1,left,bottom}
	ModifyGraph lStyle(circle_y)=1
	ModifyGraph rgb=(0,0,0)
	ModifyGraph nticks=0
	ModifyGraph font="Helvetica"
	ModifyGraph axThick=0
EndMacro