Retro Programming Nostalgia VIII: 1926-2026 l’equazione di Schrödinger e la struttura elettronica dell’atomo d’idrogeno |

Quest’anno ricorre l’anniversario della pubblicazione dell’articolo di Edwin Schrödinger (1887-1961) in cui viene introdotta la sua famosa equazione. Prendendo spunto da questa occasione, ho ripescato e rinnovato uno dei miei antichi progetti di programmazione in BASIC con i miei microcomputer negli anni ’80. Di nuovo il microcomputer era il mio amato Phillips MSX, di cui ho parlato in altri blog. Studiando chimica, non potevo non essere attratto dalla bellezza e dall’eleganza delle soluzioni dell’equazione di Schrödinger per l’atomo d’idrogeno. Inspirato dal libro (S. Marseglia, La Chimica col personal computer pubblicato dalla Muzzio) in cui mostrava alcuni esempi di programmi in BASIC per la chimica, decisi di imbarcarmi nell’impresa e usare l’MSX e poi l’Amiga Basic Basic per provare a riprodurre le bellissime visualizzazioni degli orbitali molecolari che vedevo nei libri di chimica universitari. Ma prima di questo vediamo di tornare a contenuto dell’articolo di Schrödinger.

1926–2026: cento anni dell’equazione di Schrödinger

Il 2026 segna un anniversario particolarmente significativo per la fisica moderna: cento anni dalla formulazione dell’equazione di Schrödinger, pubblicata nel 1926. Con quel lavoro, Erwin Schrödinger fornì per la prima volta una descrizione matematica coerente e continua dei sistemi quantistici, cambiando radicalmente il modo di concepire la materia a scala atomica. Fu proprio applicando l’equazione di Schrödinger all’atomo di idrogeno che emerse una delle più profonde rivoluzioni concettuali della fisica: l’elettrone non percorre orbite definite, ma è descritto da una funzione d’onda distribuita nello spazio. Le soluzioni di questa equazione portarono naturalmente alla nascita degli orbitali atomici, un’accurata rappresentazione probabilistica degli elettroni intorno agli atomi. Nel 1926 queste idee erano profondamente astratte e difficili da visualizzare. Oggi, a distanza di un secolo, possiamo invece calcolare e rappresentare graficamente le funzioni d’onda dell’atomo di idrogeno anche con semplici programmi, riscoprendo la bellezza, l’estetica e l’importanza di quella rivoluzione scientifica. Questo anno segna anche un’altra importante data, la nascita del primo personal computer commerciale, l’Apple II, messo in commercio nel 1976.

Questo blog nasce proprio con questo spirito: celebrare il centenario della meccanica ondulatoria tornando alle radici teoriche dell’atomo di idrogeno. Un piccolo omaggio al 1926, realizzato nel 2026, per ricordare come una singola equazione abbia cambiato per sempre la nostra visione della natura. Quest’anno ricorre anche il cinquantenario del personal computer. Per questo motivo vorrei fare questa celebrazione utilizzando lo stesso linguaggio di programmazione dei primi esperimenti di visualizzazione scientifica sui primi personal computers.

L’atomo di idrogeno: il laboratorio fondamentale della meccanica quantistica

L’atomo di idrogeno è l’atomo più semplice esistente in natura: è formato da un solo protone (l’isotopo più leggero, il prozio), che costituisce il nucleo, e da un solo elettrone legato ad esso dalla forza elettrostatica. Proprio per questa semplicità, l’idrogeno ha avuto un ruolo centrale nello sviluppo della meccanica quantistica ed è ancora oggi un punto di riferimento fondamentale per comprendere la struttura della materia. Dal punto di vista classico si potrebbe immaginare l’elettrone come una particella che orbita attorno al nucleo, in modo simile a un pianeta intorno al Sole. Tuttavia, questa immagine è profondamente scorretta. La descrizione corretta dell’atomo di idrogeno emerge solo all’interno della meccanica quantistica, dove lo stato dell’elettrone è rappresentato da una funzione d’onda, soluzione dell’equazione di Schrödinger. Questa funzione d’onda non descrive una traiettoria, ma una distribuzione di probabilità: il quadrato del suo modulo indica la probabilità di trovare l’elettrone in una certa regione dello spazio. Le regioni di spazio in cui questa probabilità è maggiore prendono il nome di orbitali atomici.

Gli orbitali dell’idrogeno sono caratterizzati da tre numeri quantici:

il numero quantico principale n, che determina l’energia e la dimensione dell’orbitale;
il numero quantico orbitale l, che ne determina la forma;
il numero quantico magnetico m, che ne determina l’orientazione nello spazio.

A seconda dei valori di questi numeri quantici si ottengono orbitali con forme molto diverse: sferiche (orbitali s), bilobate (p), a quattro lobi (d) o ancora più complesse (f). Una caratteristica fondamentale degli orbitali è la presenza di nodi, superfici o linee in cui la probabilità di trovare l’elettrone è esattamente nulla, una conseguenza puramente quantistica senza analoghi classici. La sezione che segue è un riassunto molto tecnico e condensato dell’equazione di Schrödinger nello stato stazionario. Se il lettore o la lettrice è interessata alla sua derivazione può lasciare un commento per stimolare un blog dedicato all’argomento. Altrimenti può consultare eccellenti testi universitari di chimica fisica che forniscono i fondamenti necessari.

Equazione di Schrödinger per l’atomo di idrogeno

L’equazione di Schrödinger stazionaria è: $\hat{H}\psi = E\psi$ . La sua energia totale (cinetica più potenziale) è calcolata usando l’operatore di Hamilton (Hamiltoniana):

$\hat{H} = -\frac{\hbar^2}{2m}\nabla^2 - \frac{e^2}{4\pi\varepsilon_0 r}$

Poiché il potenziale dipende solo da $r$ , è naturale usare coordinate sferiche. Il Laplaciano in coordinate sferiche è:

$\nabla^2 = \frac{1}{r^2}\frac{\partial}{\partial r} \left(r^2 \frac{\partial}{\partial r}\right) + \frac{1}{r^2\sin\theta} \frac{\partial}{\partial\theta} \left(\sin\theta \frac{\partial}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta} \frac{\partial^2}{\partial\phi^2}$

Separazione delle variabili

Cerchiamo soluzioni nella forma: $\psi(r,\theta,\phi) = R(r)Y(\theta,\phi)$

Ovvero come prodotto di due equazioni separate che dipendono dalle coordinate angolari (parte angolare)

$\hat{L}^2 Y = \hbar^2 l(l+1) Y$

Le soluzioni sono le armoniche sferiche:

$Y_l^m(\theta,\phi)N_{lm} P_l^m(\cos\theta) e^{im\phi}$

dove $P_l^m$ sono i polinomi di Legendre associati, l = 0,1,2,… and, m = -l,…,+l

L’altra equazione dipende dalla distanza (equazione radiale)

$\frac{d}{dr} \left(r^2 \frac{dR}{dr}\right) + \left[ \frac{2m}{\hbar^2} \left( E + \frac{e^2}{4\pi\varepsilon_0 r} \right)r^2 l(l+1) \right]R = 0$

Che viene normalizzata usando la sostituzione: $\rho = \frac{2Zr}{na_0}$ , a questa forma

$R_{nl}(r) N_{nl} e^{-\rho/2} \rho^l L_{n-l-1}^{2l+1}(\rho)$

dove $L_{n-l-1}^{2l+1}$ sono polinomi di Laguerre associati, $a_0$ è il raggio di Bohr, n = 1,2,3,…, l = 0,1,…,n-1.

Pertanto, le autofunzioni dell’atomo di idrogeno sono:

$\psi_{nlm}(r,\theta,\phi) R_{nl}(r) Y_l^m(\theta,\phi)$

Esempi di orbitali atomici

Le funzioni d’onda dell’atomo di idrogeno sono normalizzate secondo
$\int |\psi_{nlm}|^2 dV = 1$

e possono essere scritte come
$\psi_{nlm}(r,\theta,\phi) = R_{nl}(r)Y_l^m(\theta,\phi)$

dove $a_0$ è il raggio di Bohr.

1s: $\psi_{100}= \frac{1}{\sqrt{\pi a_0^3}} e^{-r/a_0}$

$2p_z: \psi_{2p_z} = \frac{1}{4\sqrt{2\pi a_0^3}} \left(\frac{r}{a_0}\right) e^{-r/2a_0} \cos\theta$

in coordinate cartesiane: $\psi_{2p_z} = \frac{1}{4\sqrt{2\pi a_0^5}} z\,e^{-r/2a_0}$

$2p_x: \psi_{2p_x} = \frac{1}{4\sqrt{2\pi a_0^5}} x\,e^{-r/2a_0}$

$2p_y: \psi_{2p_y} = \frac{1}{4\sqrt{2\pi a_0^5}} y\,e^{-r/2a_0}$

$3d_{z^2}: \psi_{3d_{z^2}} = \frac{1}{81\sqrt{6\pi a_0^7}} (3z^2 - r^2) e^{-r/3a_0}$

$3d_{x^2-y^2}: \psi_{3d_{x^2-y^2}} = \frac{1}{81\sqrt{6\pi a_0^7}} (x^2 - y^2) e^{-r/3a_0}$

$3d_{xy}: \psi_{3d_{xy}} = \frac{1}{81\sqrt{6\pi a_0^7}} xy\,e^{-r/3a_0}$

$3d_{xz}: \psi_{3d_{xz}} = \frac{1}{81\sqrt{6\pi a_0^7}} xz\,e^{-r/3a_0}$

$3d_{yz}: \psi_{3d_{yz}} = \frac{1}{81\sqrt{6\pi a_0^7}} yz\,e^{-r/3a_0}$

Visualizzare gli orbitali: dal formalismo matematico alle immagini

Come accennato nell’introduzione, il programma QB64 Basic che presento in questo blog è l’evoluzione di un programma in Basic per microcomputer che sviluppai durante gli anni universitari. Il programma utilizza le espressioni analitiche delle funzioni d’onda dell’atomo di idrogeno per calcolare le componenti radiale e angolare degli orbitali e ne rappresenta graficamente la distribuzione nello spazio.

Le immagini prodotte mostrano:

La funzione radiale (in alto a sinistra).
la funzione angolare proiettata sul piano XZ (in basso a sinistra);
la proiezione sul piano XY (in alto a destra);
la rappresentazione tridimensionale dell’isosuperficie dell’orbitale (in basso a destra).

Queste rappresentazioni permettono di evidenziare

la forma caratteristica di ciascun orbitale,
la presenza dei nodi,
l’orientazione nello spazio degli stati quantici.

Questo tipo di visualizzazione ha anche un valore storico e didattico. È lo stesso approccio che, a partire dagli anni Venti del Novecento, ha contribuito a rendere più intuitivo il concetto di densità elettronica, permettendo di collegare la teoria quantistica agli spettri atomici e alle proprietà chimiche degli elementi, e di comprendere meglio la struttura e la reattività di atomi e molecole.

Personalmente, la creazione di questo programma è stata un progetto molto piacevole, che mi ha aiutato a comprendere e demistificare la matematica apparentemente complessa dell’equazione di Schrödinger. Invito quindi tutti i lettori a esplorare e migliorare il programma in Appendice, magari riscrivendolo nel linguaggio di programmazione con cui hanno maggiore familiarità.

Appendice – Il programma

Seguendo una tradizione che era comune nelle riviste di microcomputer degli anni Settanta e Ottanta, riporto qui il listato completo del programma utilizzato per generare le visualizzazioni degli orbitali dell’atomo di idrogeno presentate in questo articolo. In quel periodo, riviste dedicate all’informatica personale pubblicavano frequentemente programmi completi che i lettori potevano digitare direttamente sul proprio computer per esplorarne il funzionamento e sperimentare modifiche.

Il programma presentato qui è scritto in QB64 Basic ed è l’evoluzione di un codice che avevo sviluppato durante gli anni universitari per visualizzare le soluzioni analitiche dell’equazione di Schrödinger per l’atomo di idrogeno. Il codice utilizza le espressioni esplicite delle funzioni radiali e degli armonici sferici per calcolare la distribuzione spaziale della funzione d’onda e rappresentare graficamente le caratteristiche principali degli orbitali atomici.

Le immagini generate dal programma includono:

la funzione radiale $R_{n\ell}(r)$ ,
la distribuzione angolare associata agli armonici sferici,
proiezioni bidimensionali della densità elettronica,
una rappresentazione tridimensionale dell’isosuperficie della funzione d’onda.

Il codice può essere utilizzato liberamente per scopi didattici. I lettori sono incoraggiati a modificarlo, esplorarne la struttura e, se lo desiderano, riscriverlo in altri linguaggi di programmazione.

Di seguito è riportato il listato completo del programma.

Nota: Il programma è fornito a scopo didattico e senza alcuna garanzia. Sebbene la versione fornita sia stata testata, l’autore non può garantire che sia priva di errori o adatta a ogni sistema.

Rem ******************************************************
Rem HYDROGENOID ORBITALS VISUALIZER
Rem QB64 Version
Rem
Rem 1926 - 2026 : 100 Years of Wave Mechanics
Rem
Rem V1.0  1987   MSX BASIC
Rem V2.0  2025-2026   Adapted for QB64 BASIC
Rem
Rem Programmed by Danilo Roccatano
Rem (c) 1987-2026
Rem
Rem This program is provided for educational purposes.
Rem ******************************************************

Screen _NewImage(800, 750, 256)
_ScreenMove 100, 50
_Display


Cls
Color 0, 15
Dim Shared PI As Double
PI = 4 * Atn(1)
a0 = 0.529 ' Bohr radius in Angstroms
ia0 = 1 / a0


Dim angY As Double
Dim angX As Double
angX = 0
angY = 0


Do

    Cls
    Call DrawFrames
    Call WriteInfo

    Color 0, 14
    Locate 35, 55: Input "Enter n, l, m "; n, l, m

    If n < 1 Or l < 0 Or l >= n Or Abs(m) > l Then
        Locate 29, 55: Print "Invalid quantum numbers"
        Locate 30, 55: Print "Press any key..."
        Do: Loop Until InKey$ <> ""
    Else

        angY = 0
        angX = 0

        ' -------- INTERACTIVE ROTATION LOOP --------
        Do

            _Limit 60
            Cls
            Call DrawFrames
            Call WriteInfo

            ' ---- Rotation keys ----
            If _KeyDown(19712) Then angY = angY - .05 ' Left
            If _KeyDown(19200) Then angY = angY + .05 ' Right
            If _KeyDown(18432) Then angX = angX + .05 ' Up
            If _KeyDown(20480) Then angX = angX - .05 ' Down

            ' Reset view
            If _KeyDown(57) Then angX = 0: angY = 0 ' Spacebar

            ' ---- Draw ----
            Call PlotRadial(n, l)
            Call PlotAngularXY(l, m)
            Call PlotAngularXZ(l, m)
            Call PlotAngular3DMeshHidden(l, m, angY, angX)
            Call DisplayQuantumInfo(n, l, m)

            Locate 44, 3: Print "Arrows = Rotate   Space = Reset   N = New orbital   ESC = Quit"

            If _KeyDown(78) Then Exit Do ' N key for new orbital
            If _KeyDown(27) Then End ' ESC quits program
            _Display
        Loop

    End If

Loop


'===========================================================
Sub DrawFrames
    Line (10, 10)-(390, 340), 15, BF
    Line (400, 10)-(780, 340), 15, BF
    Line (10, 350)-(390, 690), 15, BF
    'Line (400, 350)-(780, 690), 8, BF
    Line (405, 355)-(775, 685), 7, BF

    Line (15, 15)-(385, 335), 0, B
    Line (405, 15)-(775, 335), 0, B
    Line (15, 355)-(385, 685), 0, B
    Line (405, 355)-(775, 685), 0, B

    Color 0, 15
    'Locate 1, 4: Print "Radial probability density"
    Locate 1, 53: Print "Angular density  XY plane"
    Locate 22, 4: Print "Angular density  XZ plane"

    ' Draw axes
    Color 2, 15
    Locate 22, 45: Print "r"
    ' Add scale indicators
    Color 1, 15
    Locate 1, 4: Print "P(r) = 4pi*r^2*|R(r)|^2"
    'Locate 4, 5: Print "r in units of a0"
End Sub

'===========================================================
Sub WriteInfo
    Color 1, 7
    Locate 24, 55: Print "Hydrogen-like orbitals"
    Locate 25, 55: Print "Schrodinger equation (1926)"
    Locate 26, 55: Print "---------------------------"
    Locate 27, 55: Print "Radial:  4*pi*r^2 |R(r)|^2"
    Locate 28, 55: Print "Angular: |Y_lm|^2"
    Locate 30, 55: Print "Examples:"
    Locate 31, 55: Print " 1 0 0  -> 1s"
    Locate 32, 55: Print " 2 1 0  -> 2pz"
    Locate 33, 55: Print " 2 1 1  -> 2px"
    Locate 34, 55: Print " 2 1 -1 -> 2py"
    Locate 35, 55: Print " 3 2 0  -> 3dz^2"
End Sub

'===========================================================
Sub DisplayQuantumInfo (n, l, m)
    Color 4, 15
    Locate 38, 53: Print "Current orbital:"
    Locate 39, 53: Print "n ="; n; "  l ="; l; "  m ="; m

    ' Orbital name
    orbitalName$ = ""
    If l = 0 Then orbitalName$ = "s"
    If l = 1 Then orbitalName$ = "p"
    If l = 2 Then orbitalName$ = "d"
    If l = 3 Then orbitalName$ = "f"

    mName$ = ""
    If l = 0 Then
        mName$ = ""
    ElseIf l = 1 Then
        If m = 0 Then mName$ = "z"
        If m = 1 Then mName$ = "x"
        If m = -1 Then mName$ = "y"
    ElseIf l = 2 Then
        If m = 0 Then mName$ = "z^2"
        If m = 1 Then mName$ = "xz"
        If m = -1 Then mName$ = "yz"
        If m = 2 Then mName$ = "x^2-y^2"
        If m = -2 Then mName$ = "xy"
    End If

    Locate 40, 53: Print "Orbital: "; n; orbitalName$; mName$

    ' Radial nodes
    radialNodes = n - l - 1
    Locate 41, 53: Print "Radial nodes: "; radialNodes

    ' Angular nodes
    angularNodes = l
    Locate 42, 53: Print "Angular nodes: "; angularNodes
End Sub

'===========================================================
Sub PlotRadial (n, l)
    x0 = 25
    y0 = 320
    dr = .25

    xscale = 2.5
    yscale = 80

    ' Clear previous plot
    Line (15, 335)-(385, 335), 15, BF
    Line (15, 15)-(385, 335), 0, B

    ' Draw axes
    Line (x0, 35)-(x0, y0), 2
    Line (x0, y0)-(385, y0), 2

    ' Axis labels
    Color 1, 15
    Locate 20, 3: Print "0"
    Locate 20, 20: Print "5a0"
    Locate 20, 35: Print "10a0"

    ' Draw tick marks
    For i = 0 To 10 Step 5
        Line (x0 + i * xscale / dr, y0)-(x0 + i * xscale / dr, y0 + 5), 2
    Next i

    ' Calculate and plot radial probability density
    maxP = 0
    rValueAtMax = 0

    ' Pre-calculate to find maximum for scaling
    r = dr
    Do While r < 30
        Rnl = RadialFunction(n, l, r)

        P = 4 * PI * r * r * Rnl * Rnl

        If P > maxP Then
            maxP = P
            rValueAtMax = r
        End If
        r = r + dr
    Loop

    If maxP > 0 Then
        yscale = 250 / maxP ' Auto-scale to fit 250 pixels
    End If

    ' Now plot with proper scaling
    r = dr
    firstPoint = 1
    Do While r < 30
        Rnl = RadialFunction(n, l, r)
        P = 4 * PI * r * r * Rnl * Rnl

        xplot = x0 + xscale * (r / dr)
        yplot = y0 - yscale * P
        '     Print xplot, yplot, P, yscale
        ' Limit y to plot area
        If yplot > y0 Then yplot = y0

        If firstPoint = 1 Then
            PSet (xplot, yplot), 3
            firstPoint = 0
        Else
            If (xplot < 385) Then
                Line -(xplot, yplot), 3
            End If
        End If

        r = r + dr
    Loop

    ' Display max value
    Color 4, 15
    Locate 2, 5: Print "Max P(r) =";
    Locate 2, 20: Print Using "##.###"; maxP;
    Locate 2, 30: Print "at r =";
    Locate 2, 37: Print Using "##.##"; rValueAtMax;
    Locate 2, 42: Print "a0"
End Sub

'===========================================================
Function RadialFunction (n, l, r)
    ' Hydrogen radial wavefunctions R_nl(r)
    ' Atomic units: a0 = 1
    ' Normalized so that ∫ r^2 |R|^2 dr = 1

    Select Case n

        '-------------------------
        Case 1 ' 1s
            If l = 0 Then
                RadialFunction = 2 * Exp(-r)
            Else
                RadialFunction = 0
            End If

            '-------------------------
        Case 2 ' 2s, 2p
            If l = 0 Then
                RadialFunction = (1 / Sqr(2)) * (1 - r / 2) * Exp(-r / 2)
            ElseIf l = 1 Then
                RadialFunction = (1 / (2 * Sqr(6))) * r * Exp(-r / 2)
            Else
                RadialFunction = 0
            End If

            '-------------------------
        Case 3 ' 3s, 3p, 3d
            If l = 0 Then
                RadialFunction = (2 / (81 * Sqr(3))) * _
                                 (27 - 18 * r + 2 * r * r) * Exp(-r / 3)
            ElseIf l = 1 Then
                RadialFunction = (4 / (81 * Sqr(6))) * _
                                 r * (6 - r) * Exp(-r / 3)
            ElseIf l = 2 Then
                RadialFunction = (4 / (81 * Sqr(30))) * _
                                 r * r * Exp(-r / 3)
            Else
                RadialFunction = 0
            End If

            '-------------------------
        Case 4 ' 4s, 4p, 4d, 4f
            If l = 0 Then ' 4s
                RadialFunction = (1 / (768 * Sqr(2))) * _
                                 (192 - 144 * r + 24 * r ^ 2 - r ^ 3) * Exp(-r / 4)

            ElseIf l = 1 Then ' 4p
                RadialFunction = (1 / (768 * Sqr(6))) * _
                                 r * (80 - 20 * r + r ^ 2) * Exp(-r / 4)

            ElseIf l = 2 Then ' 4d
                RadialFunction = (1 / (768 * Sqr(30))) * _
                                 r ^ 2 * (12 - r) * Exp(-r / 4)

            ElseIf l = 3 Then ' 4f
                RadialFunction = (1 / (768 * Sqr(210))) * _
                                 r ^ 3 * Exp(-r / 4)

            Else
                RadialFunction = 0
            End If

            '-------------------------
        Case Else
            RadialFunction = 0

    End Select
End Function



'===========================================================
Sub PlotAngularXY (l, m)

    cx = 590
    cy = 170
    R0 = 120

    ' Clear previous plot
    Line (405, 15)-(775, 335), 15, BF
    Line (405, 15)-(775, 335), 0, B

    ' Draw axes
    Line (cx - R0, cy)-(cx + R0, cy), 4
    Line (cx, cy - R0)-(cx, cy + R0), 4

    ' Axis labels
    Color 4, 15
    Locate 3, 74: Print "Y"
    Locate 11, 93: Print "X"

    theta = 1.57079632679 ' π/2 → XY plane

    ' --- find max |Ylm|^2 for scaling
    maxF2 = 0
    For phi = 0 To 6.28318530718 Step 0.034906585
        F = Ylm(l, m, theta, phi)
        If F * F > maxF2 Then maxF2 = F * F
    Next phi

    If maxF2 = 0 Then Exit Sub
    ANGSCALE = R0 * 0.9 / maxF2

    ' --- draw multiple isolines
    For iso = 0.25 To 1 Step 0.25
        firstPoint = 1

        For phi = 0 To 6.28318530718 Step 0.034906585
            F = Ylm(l, m, theta, phi)
            r = ANGSCALE * iso * (F * F)

            ' +r branch
            x1 = cx + r * Cos(phi)
            y1 = cy + r * Sin(phi)

            ' −r branch (critical!)
            x2 = cx - r * Cos(phi)
            y2 = cy - r * Sin(phi)

            If F >= 0 Then
                plotColor = 2 ' green
            Else
                plotColor = 4 ' red
            End If

            If firstPoint = 1 Then
                PSet (x1, y1), plotColor
                PSet (x2, y2), plotColor
                firstPoint = 0
            Else
                Line -(x1, y1), plotColor
                Line -(x2, y2), plotColor
            End If
        Next phi
    Next iso

    ' Info text
    Color 1, 15
    Locate 19, 53: Print "XY plane (theta = 90)"
    Locate 20, 53: Print "|Y_lm|^2 isolines"

End Sub

'===========================================================

Sub PlotAngularXZ (l, m)

    cx = 200
    cy = 520
    R0 = 120

    '-------------------------------------------------
    ' Clear plot area
    '-------------------------------------------------
    Line (15, 355)-(385, 685), 15, BF
    Line (15, 355)-(385, 685), 0, B

    '-------------------------------------------------
    ' Draw axes
    '-------------------------------------------------
    Line (cx - R0, cy)-(cx + R0, cy), 4 ' X
    Line (cx, cy - R0)-(cx, cy + R0), 4 ' Z

    Color 4, 15
    Locate 24, 25: Print "Z"
    Locate 33, 45: Print "X"

    phi = 0 ' XZ plane

    '-------------------------------------------------
    ' Find max |Y_lm| for scaling
    '-------------------------------------------------
    maxF = 0
    For theta = 0 To 6.28318530718 Step 0.01 ' 0 → 2π
        F = Abs(Ylm(l, m, theta, phi))
        If F > maxF Then maxF = F
    Next theta

    If maxF = 0 Then Exit Sub
    ANGSCALE = R0 * 0.9 / maxF

    '-------------------------------------------------
    ' Draw isolines of |Y_lm|^2
    '-------------------------------------------------
    levels = 4

    For level = 1 To levels
        iso = level / levels
        firstPoint = 1

        For theta = 0 To 6.28318530718 Step 0.01

            F = Ylm(l, m, theta, phi)

            ' Signed radius → full plane
            r = ANGSCALE * iso * F * Abs(F)

            x = cx + r * Sin(theta)
            y = cy - r * Cos(theta)

            If F >= 0 Then
                plotColor = 2 + level ' green shades
            Else
                plotColor = 4 + level ' red shades
            End If

            If firstPoint = 1 Then
                PSet (x, y), plotColor
                firstPoint = 0
            Else
                Line -(x, y), plotColor
            End If

        Next theta
    Next level

    '-------------------------------------------------
    ' Info
    '-------------------------------------------------
    Color 1, 15
    Locate 41, 4: Print "XZ plane (phi = 0)"
    '    Locate 41, 4: Print "Isolines of |Y_lm|^2"
    Locate 42, 4: Print "Angular nodes:"; l

End Sub


'===========================================================
Function Ylm (l, m, theta, phi)
    ' Real spherical harmonics for hydrogen orbitals
    ' Normalized real combinations (physics / chemistry convention)
    ' theta: polar angle [0,pi]
    ' phi: azimuthal angle [0,2pi]

    Dim F As Double
    PI = 3.14159265358979

    Select Case l

        '========================
        Case 0 ' s orbital
            F = Sqr(1 / (4 * PI))

            '========================
        Case 1 ' p orbitals
            Select Case m
                Case 0 ' pz
                    F = Sqr(3 / (4 * PI)) * Cos(theta)
                Case 1 ' px
                    F = Sqr(3 / (4 * PI)) * Sin(theta) * Cos(phi)
                Case -1 ' py
                    F = Sqr(3 / (4 * PI)) * Sin(theta) * Sin(phi)
                Case Else
                    F = 0
            End Select

            '========================
        Case 2 ' d orbitals
            Select Case m
                Case 0 ' dz^2
                    F = Sqr(5 / (16 * PI)) * (3 * Cos(theta) ^ 2 - 1)

                Case 1 ' dxz
                    F = Sqr(15 / (4 * PI)) * Sin(theta) * Cos(theta) * Cos(phi)

                Case -1 ' dyz
                    F = Sqr(15 / (4 * PI)) * Sin(theta) * Cos(theta) * Sin(phi)

                Case 2 ' dx^2 - y^2
                    F = Sqr(15 / (16 * PI)) * Sin(theta) ^ 2 * Cos(2 * phi)

                Case -2 ' dxy
                    F = Sqr(15 / (16 * PI)) * Sin(theta) ^ 2 * Sin(2 * phi)

                Case Else
                    F = 0
            End Select

            '========================
        Case 3 ' f orbitals
            Select Case m
                Case 0 ' fz^3
                    F = Sqr(7 / (16 * PI)) * _
                        (5 * Cos(theta) ^ 3 - 3 * Cos(theta))

                Case 1 ' fxz^2
                    F = Sqr(21 / (64 * PI)) * _
                        Sin(theta) * (5 * Cos(theta) ^ 2 - 1) * Cos(phi)

                Case -1 ' fyz^2
                    F = Sqr(21 / (64 * PI)) * _
                        Sin(theta) * (5 * Cos(theta) ^ 2 - 1) * Sin(phi)

                Case 2 ' fz(x^2 - y^2)
                    F = Sqr(105 / (16 * PI)) * _
                        Sin(theta) ^ 2 * Cos(theta) * Cos(2 * phi)

                Case -2 ' fxyz
                    F = Sqr(105 / (16 * PI)) * _
                        Sin(theta) ^ 2 * Cos(theta) * Sin(2 * phi)

                Case 3 ' fx(x^2 - 3y^2)
                    F = Sqr(35 / (32 * PI)) * _
                        Sin(theta) ^ 3 * Cos(3 * phi)

                Case -3 ' fy(3x^2 - y^2)
                    F = Sqr(35 / (32 * PI)) * _
                        Sin(theta) ^ 3 * Sin(3 * phi)

                Case Else
                    F = 0
            End Select

            '========================
        Case Else
            F = 0

    End Select

    Ylm = F
End Function

Sub PlotAngular3DMeshHidden (l, m, angY, angX)

    panelLeft = 405
    panelRight = panelLeft + 370
    panelTop = 355
    panelBottom = 685

    panelWidth = panelRight - panelLeft
    panelHeight = panelBottom - panelTop

    cx = (panelLeft + panelRight) / 2
    cy = (panelTop + panelBottom) / 2

    alpha = 0
    beta = 0

    dTheta = .17453293
    dPhi = .08726646

    ' View direction (projection direction)
    Vx = -alpha
    Vy = -beta
    Vz = 1

    Vx = 0
    Vy = 0
    Vz = 1


    ' --- First pass: compute bounding box ---
    minX = 1E+30: maxX = -1E+30
    minY = 1E+30: maxY = -1E+30

    For theta = 0 To 3.14159265359 Step dTheta
        For phi = 0 To 6.28318530718 Step dPhi

            Call SurfacePoint(l, m, theta, phi, 350, x, y, z)

            ' Rotate Y then X
            Call RotateY(x, y, z, angY, x1, y1, z1)
            Call RotateX(x1, y1, z1, angX, xr, yr, zr)

            xs = xr + alpha * zr
            ys = yr - beta * zr

            If xs < minX Then minX = xs
            If xs > maxX Then maxX = xs
            If ys < minY Then minY = ys
            If ys > maxY Then maxY = ys


        Next phi
    Next theta

    scaleX = panelWidth / (maxX - minX)
    scaleY = panelHeight / (maxY - minY)
    scaleX = panelWidth / (maxX - minX)
    scaleY = panelHeight / (maxY - minY)

    scaleF = scaleX
    If scaleY < scaleF Then scaleF = scaleY
    scaleF = scaleF * .85
    centerX3D = (minX + maxX) / 2
    centerY3D = (minY + maxY) / 2



    Line (panelLeft, panelTop)-(panelRight, panelBottom), 15, BF
    Line (panelLeft, panelTop)-(panelRight, panelBottom), 0, B

    ' --- Second pass: draw with hidden removal ---
    For theta = 0 To 3.14159265359 - dTheta Step dTheta
        For phi = 0 To 6.28318530718 - dPhi Step dPhi

            ' --- Compute 3 surface points ---
            Call SurfacePoint(l, m, theta, phi, 350, x1, y1, z1)
            Call SurfacePoint(l, m, theta + dTheta, phi, 350, x2, y2, z2)
            Call SurfacePoint(l, m, theta, phi + dPhi, 350, x3, y3, z3)

            ' --- Rotate them ---
            Call RotateY(x1, y1, z1, angY, x1, y1, z1)
            Call RotateX(x1, y1, z1, angX, x1, y1, z1)

            Call RotateY(x2, y2, z2, angY, x2, y2, z2)
            Call RotateX(x2, y2, z2, angX, x2, y2, z2)

            Call RotateY(x3, y3, z3, angY, x3, y3, z3)
            Call RotateX(x3, y3, z3, angX, x3, y3, z3)

            ' --- Compute tangents ---
            Tx1 = x2 - x1
            Ty1 = y2 - y1
            Tz1 = z2 - z1

            Tx2 = x3 - x1
            Ty2 = y3 - y1
            Tz2 = z3 - z1

            ' --- Normal ---
            Nx = Ty1 * Tz2 - Tz1 * Ty2
            Ny = Tz1 * Tx2 - Tx1 * Tz2
            Nz = Tx1 * Ty2 - Ty1 * Tx2

            visibility = Nx * Vx + Ny * Vy + Nz * Vz

            If visibility < 0 Then

                ' --- Projection ---
                xs1 = cx + scaleF * (x1 - centerX3D)
                ys1 = cy - scaleF * (y1 - centerY3D)

                xs2 = cx + scaleF * (x2 - centerX3D)
                ys2 = cy - scaleF * (y2 - centerY3D)

                xs3 = cx + scaleF * (x3 - centerX3D)
                ys3 = cy - scaleF * (y3 - centerY3D)


                ' --- Depth fading ---
                avgZ = (z1 + z2 + z3) / 3

                depthColor = 2 + (avgZ + 350) / 700 * 10
                If depthColor < 2 Then depthColor = 2
                If depthColor > 15 Then depthColor = 15

                Line (xs1, ys1)-(xs2, ys2), depthColor
                Line (xs1, ys1)-(xs3, ys3), depthColor

            End If

        Next phi
    Next theta

End Sub

Sub SurfacePoint (l, m, theta, phi, R0, x, y, z)

    F = Ylm(l, m, theta, phi)
    r = R0 * F * F

    x = r * Sin(theta) * Cos(phi)
    y = r * Sin(theta) * Sin(phi)
    z = r * Cos(theta)

End Sub


Sub RotateY (x, y, z, ang, xr, yr, zr)
    xr = x * Cos(ang) + z * Sin(ang)
    yr = y
    zr = -x * Sin(ang) + z * Cos(ang)
End Sub

Sub RotateX (x, y, z, ang, xr, yr, zr)
    xr = x
    yr = y * Cos(ang) - z * Sin(ang)
    zr = y * Sin(ang) + z * Cos(ang)
End Sub

"… I seem […] only like a boy playing on the sea-shore, and diverting myself in now and then finding a smoother pebble or a prettier shell than ordinary, whilst the great ocean of truth lay all undiscovered before me". – Isaac Newton.

Retro Programming Nostalgia VIII: 1926-2026 l’equazione di Schrödinger e la struttura elettronica dell’atomo d’idrogeno

1926–2026: cento anni dell’equazione di Schrödinger

L’atomo di idrogeno: il laboratorio fondamentale della meccanica quantistica

Equazione di Schrödinger per l’atomo di idrogeno

Separazione delle variabili

Esempi di orbitali atomici

Visualizzare gli orbitali: dal formalismo matematico alle immagini

Leave a comment Cancel reply

1926–2026: cento anni dell’equazione di Schrödinger

L’atomo di idrogeno: il laboratorio fondamentale della meccanica quantistica

Equazione di Schrödinger per l’atomo di idrogeno

Separazione delle variabili

Esempi di orbitali atomici

Visualizzare gli orbitali: dal formalismo matematico alle immagini

Share this:

Related

Leave a comment Cancel reply