3D Representation of H Atom Orbitals

The following code generates different representations of the real H atom orbitals. The wave function and the wave function squared are represented as isosurfaces in different Gizmo plots. The sign of the wave function is represented by color (i.e., blue and red isosurfaces).  In a different window, a Monte Carlo simulation is used to generate a probabilistic point cloud representation of the wave function. Graphs of the radial wave function are also generated.

The isosurfaces may need to be tweaked. For people not used to the Gizmo interface, this can be accomplished with the Modify Isosurface macro.

Code to generate the Bohr orbitals is also included.

Sample images are shown after the code.

#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