3D Representation of H Atom Orbitals

#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