! 水素原子（1電子原子）の原子軌道の表示プログラム
!
!
! （十進BASIC ソース）：　全文コピーして実行してください。
!
!   原子軌道の形などに関して、図をウェブや書籍に頼るのではなく、
!   好きな時に好きな形式で表示することで、学習者のより深い理解につながると考え、
!   十進ＢＡＳＩＣを用いて、表示用のプログラムを記述してみました。
! 
!   煩雑な原子軌道についての数式は、プログラム中に定義しておきましたので、
!   比較的単純な指定だけで原子軌道を様々な方法で二次元に表示します。
!   また、自分で軌道を混成することもできます（厳密ではありません）。
!
!   商用でなければ、自由に利用したり、解析してみてください。
!   サブルーチンは、すべて内部副プログラムとして記述しています。
!   そのため、局所的な媒介変数を除いて、変数は無条件に引き継がれます。
!   動作確認環境は、十進BASIC の ver 7.83 / Win 8.1 Pro です。
!   ※ 初期に設定されているグラフィックス画面のアスペクト比を知るために
!      将来廃止予定の [ASK BITMAP SIZE x,y] を使用しています。 
!
!   十進BASIC は、以下のURIより入手してください。
!   http://hp.vector.co.jp/authors/VA008683/
!  
! （註）原子軌道は、z軸（表示上のy軸）に対する回転体として定義されています。
!       そのため、表示上のy軸に左右対称ではない混成軌道では、
!       等高線表示では、厳密な数値に基づく線ではなくなります。
!   
!
!                                                       2019.11.27  Opened
!                                                       2022.01.26  Updated
!                                                       S. Nakajima
!
DECLARE FUNCTION WFSH2
DECLARE FUNCTION WF2
! ================================================================================================
! 原子軌道の選択のための GUI
! 一部を除く混成軌道の表示、動径分布関数の表示等については、GUI からは選択できません。
CALL GUI 

! ============================================================================
! 波動関数（動径波動関数×球面調和関数）/ 球面調和関数等の指定と表示（GUI をオフの時）
! [ 使い方 ] [■] を引用して ● の関数を定義、必要なサブルーチンを CALL する。
! CALL 文（サブルーチンの呼び出し）は、不要はものはコメントアウトして下さい。
! DEF 文はコメントアウトできません。

DEF WF(r,t) = (R4s(r)*Ts0(t) - R4p(r)*Tp0(t)*SQR(3))/2   ! ● 波動関数（動径波動関数×球面調和関数）を指定（左例は、4sp3）

! 電子の存在確率に応じた大きさの円を格子点上に描く
! CALL latice

! 電子の存在確率に応じて、ランダムに点を打つ
CALL MonteCarlo

! 電子密度に応じた等高線をひく。
! LET conturLaticeSize = 0.5                         ! 大きいと、等高線の一部に抜けが生じる。小さくすればより精緻になりますが、時間が掛かる。指定しない場合、標準 1。
CALL ContourBasedOnInteg90                           ! 電子存在確率が 90 % となる領域
! CALL ContourBasedOnInteg                           ! 電子存在確率 10 % から 90 % まで10ポイント刻み 

! ============================
! 球面調和関数の形の表示
DEF WFsh(t) = (-Ts0(t) - Tp0(t)*SQR(3))/2            ! ● 球面調和関数を指定（波動関数と独立に指定することが必要）

! 定義例
! (Ts0(t) + Tp0(t))/SQR(2)                           ! sp 混成軌道-1
! (Ts0(t) - Tp0(t))/SQR(2)                           ! sp -2
! (Ts0(t) + Tp0(t)*SQR(2))/SQR(3)                    ! sp2-1
! (Ts0(t)*SQR(2) - Tp0(t) + Tp1(t)*SQR(3))/SQR(6)    ! sp2-2
! (Ts0(t)*SQR(2) - Tp0(t) - Tp1(t)*SQR(3))/SQR(6)    ! sp2-2
! (Ts0(t) + Tp0(t)*SQR(3))/SQR(4)                    ! sp3-1
! (Ts0(t)*SQR(3) - Tp0(t) + Tp1(t)*sqr(8) )/SQR(12)  ! sp3-2
! (Ts0(t)*SQR(3) - Tp0(t) - Tp1(t)*sqr(2) )/SQR(12)  ! sp3-3,4（面外成分は式に含まれない）

! 球面調和関数の二乗（*4π）を偏角 tに対してプロットする
! CALL SHSqShape2

! 球面調和関数の二乗を、次の絶対値（※）に合わせたスケールで偏角 tに対してプロットする
! CALL SHSqShape

! 球面調和関数（*2√π) の絶対値（※）を偏角 tに対してプロットする
! CALL SHShape

! 円周上に球面調和関数の絶対値に比例する径の円を並べる
! CALL SHonSphere

! ============================
! 動径波動関数、動径分布関数の表示（関数定義不要）

! 動径波動関数 wf を表示する
! CALL PlotRWF
! CALL PlotRWF_K
! CALL PlotRWF_L
! CALL PlotRWF_M
! CALL PlotRWF_N

! 動径分布関数 wf^2×r^2 を表示する
! CALL PlotRDF
! CALL PlotRDF_s
! CALL PlotRDF_p
! CALL PlotRDF_d
! CALL PlotRDF_f
! CALL plotRDF_K
! CALL plotRDF_L
! CALL plotRDF_M
! CALL plotRDF_N

! ============================================================================
! 関数定義と配列の宣言 [■] ： ●の定義用
!   ただし、表示するのは、zx 平面での断面図。
!   なお、プログラム内では直交座標に直す際、z 軸を y で表現している。

! 動径波動関数定義： 添え字の数字は主量子数、spdf は方位量子数より決まる軌道の形
DEF Z  = 1    ! 核電荷                   ' let 文として書く場合は、call より前に置く必要あり
DEF A0 = 1    ! ボーア半径を 1 とする。  ' let 文として書く場合は、call より前に置く必要あり
DEF rho(r) = Z*r/a0
DEF R1s(r) = (2) * (Z/a0)^(3/2)  * EXP(-rho(r))
DEF R2s(r) = (1/2/SQR(2)) * (Z/a0)^(3/2) * (2 - rho(r)) * EXP(-rho(r)/2)
DEF R3s(r) = (2/81/SQR(3)) * (Z/a0)^(3/2) * (27 - 18*rho(r) + 2*rho(r)^2) * EXP(-rho(r)/3)
DEF R4s(r) = (1/768) * (Z/a0)^(3/2) * (192 - 144*rho(r) + 24*rho(r)^2 - rho(r)^3) * EXP(-rho(r)/4)
DEF R2p(r) = (1/2/SQR(6)) * (Z/a0)^(3/2) * rho(r) * EXP(-rho(r)/2)
DEF R3p(r) = (4/81/SQR(6)) * (Z/a0)^(3/2) * (6 - rho(r))*rho(r) * EXP(-rho(r)/3)
DEF R4p(r) = (1/256/SQR(15)) * (Z/a0)^(3/2) * (80 - 20*rho(r) + rho(r)^2)*rho(r) * EXP(-rho(r)/4)
DEF R3d(r) = (4/81/SQR(30)) * (Z/a0)^(3/2) * Z^2*r^2/a0^2 * EXP(-rho(r)/3)
DEF R4d(r) = (1/768/SQR(5)) * (Z/a0)^(3/2) * (12 - rho(r))*rho(r)^2 * EXP(-rho(r)/4)
DEF R4f(r) = (1/768/SQR(35)) * (Z/a0)^(3/2) * rho(r)^3 * EXP(-rho(r)/4)

! 球面調和関数（二次元）定義：spdf の後の数字は磁気量子数
! z軸（表示上のy軸）に対する回転体で、磁気量子数 0 以外は2つの縮退した軌道です。
DEF Ts0(t)  = 1/2/SQR(PI)
DEF Tp0(t) = SQR(3/PI)/2 * SIN(t)
DEF Tp1(t) = SQR(3/PI)/2 * COS(t)
DEF Td0(t) = SQR(5/PI)/4 * (3*SIN(t)^2 - 1)
DEF Td1(t) = SQR(15/PI)/2 * SIN(t)*COS(t)
DEF Td2(t) = SQR(15/PI)/4 * COS(t)^2
DEF Tf0(t) = SQR(7/PI)/4 * (5*SIN(t)^3 - 3*SIN(t))
DEF Tf1(t) = SQR(10.5/PI)/4 * (5*SIN(t)^2 - 1)*COS(t)
DEF Tf2(t) = SQR(105/PI)/4 * SIN(t)*COS(t)^2
DEF Tf3(t) = SQR(17.5/PI)/4 * COS(t)^3

! 等高線表示用のパラメータ格納用の配列変数
DIM LX(4)     ! 格子の頂点のx座標
DIM LY(4)     ! 格子の頂点のy座標
DIM LH(4)     ! 各頂点における波動関数の二乗の値
DIM LS(4)     ! 各頂点における波動関数の符号

! ============================================================================
! 以下、サブルーチン

! フラグなどの初期化
sub init
   let FlagSH = 0
   let PlFlag = 0
   LET WFSM = 0  
end sub

! 電子密度、等高線等の表示用
SUB window1
   ASK bitmap SIZE a,b           ! ASK PIXEL SIZE (0 , 0 ; 1 , 1) a,b
   LET corf = (a+1)/(b+1)        ! LET corf = a/b   , 画面表示のアスペクト比
   IF w = 0 THEN LET w = 40      ! ボーア半径を 1 とする単位で定める表示範囲限界、計算範囲限界にも流用
   SET axis COLOR 1
   SET WINDOW -w*corf,w*corf,-w,w
   IF w >10 THEN
      DRAW axes (10,10)
   ELSE
      DRAW axes (2,2)
   END IF
END sub

! 球面調和関数の表示用
SUB window2
   IF w2 = 0 THEN LET w2 = 3.5   ! 
   ASK bitmap SIZE a,b
   LET corf = (a+1)/(b+1)
   SET axis COLOR 1
   SET WINDOW -w2*corf, w2*corf, -w2, w2
END SUB

! ============================================================================
! 動径波動関数の表示
SUB PlotRWF
   LET w = 40
   SET WINDOW -1, w, -0.12, 0.6
   SET axis COLOR 1
   DRAW axes (5,0.2) 
    
   SET LINE COLOR 3
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4f(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 4
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3d(x);
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4d(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 2
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2p(x);
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3p(x);
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4p(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 1
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R1s(x);
   NEXT x
   PLOT lines
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2s(x);
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3s(x);
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4s(x);
   NEXT x
   PLOT LINES
    
END SUB

SUB PlotRWF_K
   LET w = 40
   SET WINDOW -1, w, -0.4, 2
   SET axis COLOR 1
   DRAW axes (5,0.6) 
    
   SET LINE COLOR 1
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R1s(x);
   NEXT x
   PLOT lines
    
END SUB

SUB PlotRWF_L
   LET w = 40
   SET WINDOW -1, w, -0.12, 0.6
   SET axis COLOR 1
   DRAW axes (5,0.2) 
    
   SET LINE COLOR 2
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2p(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 1
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2s(x);
   NEXT x
   PLOT lines
    
END SUB

SUB PlotRWF_M
   LET w = 40
   SET WINDOW -1, w, -0.06, 0.3
   SET axis COLOR 1
   DRAW axes (5,0.1) 
    
   SET LINE COLOR 4
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3d(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 2
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3p(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 1
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3s(x);
   NEXT x
   PLOT lines
    
END SUB

SUB PlotRWF_N
   LET w = 40
   SET WINDOW -1, w, -0.05, 0.25
   SET axis COLOR 1
   DRAW axes (5,0.1) 
    
   SET LINE COLOR 3
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4f(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 4
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4d(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 2
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4p(x);
   NEXT x
   PLOT LINES
    
   SET LINE COLOR 1
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4s(x);
   NEXT x
   PLOT lines
    
END SUB

! 動径分布関数の表示
SUB PlotRDF
   LET w = 40
   SET WINDOW -1, w, -0.06, 0.6
   SET axis COLOR 1
   DRAW axes (5,0.2) 
   LET PlFlag = 1
    
   SET LINE COLOR 3
   CALL plotRDF_f 
    
   SET LINE COLOR 4
   CALL plotRDF_d
    
   SET LINE COLOR 2
   CALL plotRDF_p
    
   SET LINE COLOR 1
   CALL plotRDF_s        
END SUB

SUB plotRDF_s
   LET w = 40
   SET WINDOW -1, w, -0.06, 0.6
   SET axis COLOR 1
   DRAW axes (5,0.2) 
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R1s(x)^2*x^2;
   NEXT x
   PLOT lines
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2s(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3s(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4s(x)^2*x^2;
   NEXT x
   PLOT LINES
END sub

SUB PlotRDF_p
   IF plFlag <> 1 THEN
      LET w = 40
      SET WINDOW -1, w, -0.03, 0.3
      SET axis COLOR 1
      DRAW axes (5,0.1) 
   END IF
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2p(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3p(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4p(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_d
   IF plFlag <> 1 THEN
      LET w = 40
      SET WINDOW -1, w, -0.015, 0.15
      SET axis COLOR 1
      DRAW axes (5,0.05) 
   END IF
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3d(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4d(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_f
   IF plFlag <> 1 THEN
      LET w = 40
      SET WINDOW -1, w, -0.009, 0.09
      SET axis COLOR 1
      DRAW axes (5,0.03) 
   END IF
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4f(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_K
   LET w = 40
   SET WINDOW -1, w, -0.06, 0.6
   SET axis COLOR 1
   DRAW axes (5,0.2)
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R1s(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_L
   LET w = 40
   SET WINDOW -1, w, -0.03, 0.3
   SET axis COLOR 1
   DRAW axes (5,0.1)
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2s(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R2p(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_M
   LET w = 40
   SET WINDOW -1, w, -0.015, 0.15
   SET axis COLOR 1
   DRAW axes (5,0.05)
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3s(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3p(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R3d(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

SUB plotRDF_N
   LET w = 40
   SET WINDOW -1, w, -0.01, 0.10
   SET axis COLOR 1
   DRAW axes (5,0.1/3)
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4s(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4p(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4d(x)^2*x^2;
   NEXT x
   PLOT LINES
    
   FOR x = 0 TO w STEP 0.1
      PLOT LINES : x, R4f(x)^2*x^2;
   NEXT x
   PLOT LINES
END SUB

! ============================================================================
! 球面調和関数の表示

! 球面調和関数を動径とし、偏角に対しプロットする
SUB SHShape
   IF FlagSH<>3 THEN 
      LET w2 = 3.5
      CALL window2 
      LET FlagSH = 1
   END IF
   FOR t = 0 TO 2*PI STEP 0.01
      IF SGN(WFsh2(t))>0 THEN 
         SET LINE COLOR 2
      ELSE
         SET line COLOR 4
      END IF
      LET r = ABS(WFsh2(t)) *2*SQR(PI)
      PLOT LINES : r*COS(t), r*SIN(t);
   NEXT t
   PLOT LINES 
END SUB

! 球面調和関数の二乗を動径とし、偏角に対しプロットする（最大値を SHShape に合わせる）
SUB SHSqShape
   IF FlagSH<>3 THEN 
      LET w2 = 3.5
      CALL window2 
      LET FlagSH = 2
   END IF
   LET rmax = 0
   FOR t = 0 TO 2*PI STEP 0.01
      IF rmax < ABS(WFsh2(t)) THEN LET rmax = ABS(WFsh2(t)) 
   NEXT t
    
   FOR t = 0 TO 2*PI STEP 0.01
      IF SGN(WFsh2(t))>0 THEN 
         SET LINE COLOR 2
      ELSE
         SET line COLOR 4
      END IF
      LET r = WFsh2(t)^2/rmax *2*SQR(PI)
      PLOT LINES : r*COS(t), r*SIN(t);
   NEXT t
   PLOT LINES 
END SUB

! 球面調和関数の二乗を動径とし、偏角に対しプロットする
SUB SHSqShape2
   LET w2 = 10
   CALL window2 
   LET FlagSH = 3    
    
   FOR t = 0 TO 2*PI STEP 0.01
      IF SGN(WFsh2(t))>0 THEN 
         SET LINE COLOR 2
      ELSE
         SET line COLOR 4
      END IF
      LET r = WFsh2(t)^2*4*PI
      PLOT LINES : r*COS(t), r*SIN(t);     
   NEXT t 
   PLOT LINES 
    
END SUB

! 球面調和関数の絶対値に比例した径の円を同心円上に表示する
SUB SHonSphere
   LET w2 = 3.5
   CALL window2 
   LET rmax = 0
   FOR t = 0 TO 2*PI STEP 0.01
      IF rmax < ABS(WFsh2(t)) THEN LET rmax = ABS(WFsh2(t)) 
   NEXT t
   SET POINT STYLE 4
   FOR t = 0 TO 2*PI STEP PI/24
      LET r0 = ABS(WFsh2(t))  
      LET r1 = ABS(WFsh2(t)) *2*SQR(PI)
      LET r2 = WFsh2(t)^2/rmax *2*SQR(PI)
      LET r3 = WFsh2(t)^2*4*PI
      SELECT CASE FlagSH
      CASE 1
         LET r = r1
         LET rd = 3
      CASE 2
         LET r = r2
         LET rd = 3
      CASE 3
         LET r = r3
         LET rd = 9
      CASE ELSE
         LET r = 0
         LET rd = w2*0.9
      END SELECT
       
      IF SGN(WFsh2(t))>0 THEN 
         SET LINE COLOR 2
         SET AREA COLOR 2
      ELSE
         SET line COLOR 4
         SET AREA COLOR 4
      END IF
       
      SET LINE STYLE 1
      PLOT LINES : 0,0; COS(t)*r, SIN(t)*r 
      PLOT POINTS : COS(t)*r, SIN(t)*r ; COS(t)*r, SIN(t)*r    
      DRAW disk WITH SCALE(SQR(r0)*rd/10)*SHIFT(COS(t)*rd, SIN(t)*rd)     
      SET LINE STYLE 3
      IF FlagSH <> 0 THEN PLOT LINES : COS(t)*r, SIN(t)*r ; COS(t)*rd, SIN(t)*rd
   NEXT t 
END SUB

! ============================================================================
! 動径波動関数×球面調和関数の表示

! x,y 座標を極座標 r, t に変換する
SUB ConvertXyToRt
   LET r = SQR(x^2+y^2)
   IF x^2+y^2 = 0 THEN 
      LET t = 0
   ELSE
      LET t = SGN(y)*ACOS(x/SQR(x^2+y^2))
   END if
END SUB

! 乱数で座標を発生させ、存在確率（波動関数の二乗）と比較して点を打つ
SUB MonteCarlo
   CALL window1
   Call WFSMax
   SET POINT STYLE 1
   LET NP = 20000                ! プロットする点の数の上限、画面サイズや軌道の種類により変更可
   LET NC = 10000000             ! 試行回数の上限。
   LET c = 0  
   LET k = 0
   DO
      LET x = RND*2*w-w          ! アスペクト比が 1 ではない場合も w より外側は試行しない。
      LET y = RND*2*w-w
      CALL ConvertXyToRt         ! x,y 座標値を極座標 r,t に変換
       
      IF RND < WF2(r,t)^2/WFSM THEN
         IF SGN(WF2(r,t))>0 THEN 
            SET POINT  COLOR 2
         ELSE
            SET POINT COLOR 4
         END IF
         PLOT POINTS : x, y
         LET k = k + 1
      END IF
      LET c = c + 1
      IF k => NP THEN EXIT DO
      IF c => NC THEN EXIT DO
   loop  
   ! PRINT "試行回数"; c ; ", 打点数 ";  k
END SUB

! 格子点上に波動関数の対数に比例した半径をもつ円を描く
SUB latice
   CALL window1  
   Call WFSMax
   LET d  = 1                          ! 格子間隔
   LET df = 3                          ! 数字が大きいほど波動関数値の小さいところまでプロットする
   FOR    x = (-w*corf)+d/2 TO w*corf STEP d
      FOR y = -w+d/2 TO w STEP d
         CALL ConvertXyToRt
         IF SGN(WF2(r,t))>0 THEN 
            SET AREA COLOR 2
         ELSE
            SET AREA COLOR 4
         END if
         IF WF2(r,t) = 0 THEN            ! 計算上、log10(0) を回避するための措置。 
            LET dp = 0
         ELSE
            LET dp = MAX(0, (1 + LOG10( WF2(r,t)^2/WFSM )/df)/2*d)
         END IF
         DRAW disk WITH SCALE(dp)*SHIFT(x,y)
      NEXT y
   NEXT x
END SUB

! 波動関数の二乗値に対し、閾値 H の値の等高線を表示するための描画エンジン
SUB contour   
   LET d = 1                    ! 等高線の有無を判断するためのマス眼のサイズ
   IF conturLaticeSize  = 0 THEN               
   ELSE                         ! 外部から conturLaticeSize として指定した場合は、反映させる
      LET d = conturLaticeSize 
   END IF
   IF QNN = 1 OR w<=10 THEN LET d = 0.25     
   FOR    dx = (-w*corf)-d/2 TO w*corf STEP d
      FOR dy = -w-d/2 TO w STEP d
       
         LET LF = 1   
         FOR    x = dx TO dx+d STEP d    ! 格子の4頂点で、波動関数の二乗値と符号を調べる
            FOR y = dy TO dy+d STEP d
               CALL ConvertXyToRt
               LET LH(LF) = WF2(r,t)^2
               LET LS(LF) = SGN(WF2(r,t))           
               LET LF = LF + 1
            NEXT y
         NEXT x
          
         LET flag = 0
          
         ! x軸、y軸の周辺は細かく見る。 
         IF ABS(x) <= 2*d and ABS(y) <= 2*d   THEN LET flag = 1
          
         ! 隣接頂点間で、波動関数の二乗値と閾値の上下関係が変化する場合
         IF SGN(LH(1)-H) <> SGN(LH(2)-H) THEN LET flag = flag + 1
         IF SGN(LH(1)-H) <> SGN(LH(3)-H) THEN LET flag = flag + 1
         IF SGN(LH(2)-H) <> SGN(LH(4)-H) THEN LET flag = flag + 1
         IF SGN(LH(3)-H) <> SGN(LH(4)-H) THEN LET flag = flag + 1
          
         ! 隣接頂点間に、波動関数の節がある場合      
         IF LS(1) <> LS(2) THEN LET flag = flag + 1
         IF LS(1) <> LS(3) THEN LET flag = flag + 1
         IF LS(2) <> LS(4) THEN LET flag = flag + 1
         IF LS(3) <> LS(4) THEN LET flag = flag + 1
          
         IF flag > 1 THEN                                 ! 一辺の長さ d の格子内に線を引く        
         ! ^--------------------------------------------------------------------------------
            LET dd = d/33                                 ! 格子を更に再分割
            FOR    dx2 = dx-dd TO dx+d+dd STEP dd
               FOR dy2 = dy-dd TO dy+d+dd STEP dd
                
                  LET LF = 1   
                  FOR    x = dx2 TO dx2+dd STEP dd        ! 格子の4頂点で、波動関数の二乗値を調べる
                     FOR y = dy2 TO dy2+dd STEP dd
                        CALL ConvertXyToRt
                        LET LH(LF) = WF2(r,t)^2
                        LET LF = LF + 1
                     NEXT y
                  NEXT x
                   
                  LET flag = 0
                  IF SGN(LH(1)-H) <> SGN(LH(2)-H) THEN    ! 隣接格子間で、閾値との上下が変化したら
                     LET flag = flag + 1                  ! フラグをたて、座標 Lx, Ly を計算する。
                     LET lx(flag) = dx2
                     LET ly(flag) = dy2 + dd*LH(1)/(LH(2)+LH(1))
                  END IF
                  IF SGN(LH(1)-H) <> SGN(LH(3)-H) THEN
                     LET flag = flag + 1
                     LET lx(flag) = dx2 + dd*LH(1)/(LH(3)+LH(1))
                     LET ly(flag) = dy2
                  END IF
                  IF SGN(LH(2)-H) <> SGN(LH(4)-H) THEN
                     LET flag = flag + 1
                     LET lx(flag) = dx2 + dd*LH(2)/(LH(2)+LH(4))
                     LET ly(flag) = dy2 + dd
                  END IF
                  IF SGN(LH(3)-H) <> SGN(LH(4)-H) THEN
                     LET flag = flag + 1
                     LET lx(flag) = dx2 + dd
                     LET ly(flag) = dy2 + dd*LH(3)/(LH(3)+LH(4))
                  END IF
                   
                  LET x = dx2+dd/2
                  LET y = dy2+dd/2
                  CALL ConvertXyToRt
                   
                  IF WF2(r,t) > 0 THEN    ! 格子中央での波動関数の符号により線色を変える
                     SET LINE COLOR 2
                  ELSE 
                     SET LINE COLOR 4
                  END if
                   
                  IF flag = 2 THEN   PLOT LINES : lx(1),ly(1); lx(2),ly(2)
                   
                  IF flag = 4 THEN
                     PLOT LINES : lx(3),ly(3); lx(4),ly(4)
                     ! PRINT "【確認して下さい】 x, y ="; x; y; "で、格子内に2本の線を引きました。"
                     ! ／／ と引くべきか、＼＼ と引くべきかの二択になる。
                     ! このプログラムでは、判定せず、＼＼ と引いている。
                     ! 格子一辺の長さ dd が十分に小さければ、画面上さほど目立たない。 
                  END if  
                   
               NEXT dy2
            NEXT dx2                       
            ! ---------------------------------------------------------
         END if  
      NEXT dy
   NEXT dx
END SUB

! 波動関数の二乗の最大値の探索
SUB WFSMax
   IF WFSM = 0 THEN  ! ● 同じ軌道で、等高線を複数書くときに、再計算しなくて良いように指定したが、GUI では軌道種類を変えるごとにクリア必要
      FOR r = 0 TO 20 STEP 0.01
         FOR t = 0 TO 2*PI STEP 0.01
            IF WFSM < WF2(r,t)^2 THEN LET WFSM = WF2(r,t)^2
         NEXT t
      NEXT r
      ! PRINT "WFSM =" ; WFSM             ! 印字出力
   end if
END SUB

! 波動関数の二乗に対し等高線を引く
SUB contourBasedOnWFSM
   CALL window1 
   CALL WFSMax
   FOR k = -0.1 TO -3.1 STEP -0.5   ! WFSM（波動関数の二乗の最大値）の 10^(k) 倍を閾値とする
      LET H = WFSM * 10^(k)         ! 閾値。波動関数の二乗が H であるところに等高線を引く
      CALL contour
   NEXT k
END SUB

! 電子密度に応じた波動関数の二乗の等高線を引く。
! y軸に対して左右対称（または点対称）ではない軌道に対しては厳密ではない。
SUB contourBasedOnInteg
   CALL window1  
   Call WFSMax
   CALL INtegral                 ! y軸に対して左右対称、点対称以外の軌道では意味のある数値にならない
    
   LET partialInt = 0            ! 波動関数の二乗が大きい区間から、回転体としての微小体積を足していく
   LET partialIntApre = 0        ! 区間の数値を足す前の値
   LET dr = 0.1 
   LET dk = 0.1
   FOR k = 0 TO 10 STEP dk       ! 区間区切り位置、波動関数の二乗値の最大値に対する対数値の負数
    
      FOR r = dr/2 TO MAX(40,w) STEP dr
         FOR t = -PI/2 TO PI/2 STEP PI/180
            IF WFSM*10^(-K-dk) < WF2(r,t)^2 AND WF2(r,t)^2 <= WFSM*10^(-K) THEN
               LET partialInt = partialInt + WF2(r,t)^2*2*PI*(r*COS(t))*r*dr*PI/180
                
            END if
         NEXT t 
      NEXT r
      ! PRINT k; WFSM*10^(-K); partialIntApre ; partialInt
       
      FOR threshold = 10 TO 90 STEP 10   ! 等高線を描く電子密度の百分率
         call calcH
      NEXT threshold
      LET threshold = 90                 ! Next で 100 になっているので戻す。
       
      ! let threshold = 95               ! 等高線を描く電子密度の百分率
      ! call calcH
       
      ! let threshold = 99               ! 等高線を描く電子密度の百分率
      ! call calcH
       
      ! let threshold = 99.9             ! 等高線を描く電子密度の百分率
      ! call calcH
       
      LET partialIntApre = partialInt
      IF partialInt > roundedInt * threshold/100 THEN EXIT FOR    
   NEXT k
END sub

! 電子密度90%の等高線を1本引く
SUB contourBasedOnInteg90
   CALL window1  
   Call WFSMax  
   CALL INtegral                         ! y軸に対して左右対称、点対称以外の軌道では意味のある数値にならない
    
   LET partialInt = 0                    ! 波動関数の二乗が大きい区間から、回転体としての微小体積を足していく
   LET partialIntApre = 0                ! 区間の数値を足す前の値
   LET dr = 0.1 
   LET dk = 0.1
   FOR k = 0 TO 10 STEP dk               ! 区間区切り位置、波動関数の二乗値の最大値に対する対数値の負数
    
      FOR r = dr/2 TO MAX(40,w) STEP dr
         FOR t = -PI/2 TO PI/2 STEP PI/180
            IF WFSM*10^(-K-dk) < WF2(r,t)^2 AND WF2(r,t)^2 <= WFSM*10^(-K) THEN
               LET partialInt = partialInt + WF2(r,t)^2*2*PI*(r*COS(t))*r*dr*PI/180
                
            END if
         NEXT t 
      NEXT r
       
      LET threshold = 90                 ! 等高線を描く電子密度の百分率
      call calcH
       
      LET partialIntApre = partialInt
      IF partialInt > roundedInt * threshold/100 THEN EXIT FOR 
   NEXT k
END SUB

SUB calcH                                ! ContourBasedOnInteg 用のサブルーチン
   IF partialIntApre/roundedInt < threshold/100 AND threshold/100 < partialInt/roundedInt THEN
      LET Ktemp = k + dk*(threshold/100 - partialIntApre/roundedInt)/(partialInt/roundedInt - partialIntApre/roundedInt)
      LET H = WFSM * 10^(-Ktemp)
      PRINT "全電子密度の上位 " & STR$(threshold) & " % の等高線"  ! これからひく線の  , H/WFSM
      ! pause
      CALL contour
   END IF
End sub

! y軸に対する回転体としての微小体積の電子存在確率（Wf^2 × 2π(r×cosθ) ×dr×rdθ）の総和計算
! 全体積に対する電子存在確率の和 integ は m=0 の軌道で 1、m= ±1,2,… の軌道で 2 となる。
! ただし、y軸に対して左右対称（または点対称）ではない軌道に対しては意味のない数値を与える。
! ☆ の出力結果を、エクセルで降順に並べ、電子存在確率の積分がある特定の値に達する閾値を調べることもできる。
SUB INtegral
   LET integ = 0                             ! 微小体積ごとの値の総和
   LET dr = 0.1
   FOR r = dr/2 TO MAX(40,w) STEP dr
      FOR t = -PI/2 TO PI/2 STEP PI/180
         LET wfsitemp = WF2(r,t)^2*2*PI*(r*COS(t))*r*dr*PI/180
         ! PRINT WF2(r,t)^2; ","; wfsitemp    ! ☆ 波動関数の二乗値、微小体積ごとの値を印字（カンマ区切りテキスト）
         LET integ = integ + wfsitemp
      NEXT t 
   NEXT r
   ! PRINT integ                             ! 微小体積ごとの値の総和を印字
   LET roundedInt = ROUND(integ,2)           ! 小数点以下2桁に丸めたもの。左右非対称な軌道では小さい値をとり得る。
END SUB

! ============================================================================
! 関数定義の分岐… 条件分岐中であっても、同じ関数名で2回以上Def文で記述できないため

FUNCTION WFsh2(t)  ! FOR 
   SELECT CASE GUIflag
   CASE 1            ! GUIからの関数呼び出しである場合
      SELECT CASE QNM
      CASE 1 ! s軌道
         LET WFsh2 = Ts0(t)
      CASE 2 ! p軌道
         SELECT CASE QNL
         CASE 0
            LET WFsh2 = Tp0(t)
         CASE 1
            LET WFsh2 = Tp1(t)
         END SELECT
      CASE 3 ! d軌道
         SELECT CASE QNL
         CASE 0
            LET WFsh2 = Td0(t)
         CASE 1
            LET WFsh2 = Td1(t)
         CASE 2
            LET WFsh2 = Td2(t)
         END SELECT
      CASE 4 ! f軌道
         SELECT CASE QNL
         CASE 0
            LET WFsh2 = Tf0(t)
         CASE 1
            LET WFsh2 = Tf1(t)
         CASE 2
            LET WFsh2 = Tf2(t)
         CASE 3
            LET WFsh2 = Tf3(t)
         END SELECT
      CASE ELSE
      END SELECT
   CASE 2  ! GUI2からの呼び出し
    
      SELECT CASE FnFlag                 ! ●●●●●
       
      CASE 1                             ! ● 混成 p
         LET WFsh2 = SQR(FactorA/100)*Tp0(t) + SQR((100-FactorA)/100)*Tp1(t)
          
      CASE 2                             ! ● sp
         LET WFsh2 = (Ts0(t) + Tp0(t))/SQR(2)
          
      CASE 3                             ! ● sp2-1
         LET WFsh2 = Ts0(t)/SQR(3)  + Tp0(t)*SQR(2/3)         
          
      CASE 4                             ! ● sp2-2
         LET WFsh2 = Ts0(t)/SQR(3)  - Tp0(t)/SQR(6) + Tp1(t)/SQR(2) 
          
      CASE 5                             ! ● sp2-3
         LET WFsh2 = Ts0(t)/SQR(3)  - Tp0(t)/SQR(6) - Tp1(t)/SQR(2)
          
      CASE 6                             ! ● sp3-1
         LET WFsh2 = Ts0(t)/2  + Tp0(t)*SQR(3/4)         
          
      CASE 7                             ! ● sp3-2
         LET WFsh2 = (Ts0(t)*SQR(3) - Tp0(t) + Tp1(t)*sqr(8))/SQR(12) 
          
      CASE 8                             ! ● sp3-3
         LET WFsh2 = (Ts0(t)*SQR(3) - Tp0(t) - Tp1(t)*sqr(2))/SQR(12) 
          
      CASE 9                             ! ● sp5-1
         LET WFsh2 = (Ts0(t) + Tp0(t)*SQR(3) + Tp1(t)*sqr(2))/SQR(6) 
          
      CASE 10                            ! ● sp5-2
         LET WFsh2 = (Ts0(t) - Tp0(t)*SQR(3) + Tp1(t)*sqr(2))/SQR(6) 
          
      CASE 11                            ! ● dsp2　 3d2 4s 4p1  正方平面型
         LET WFsh2 = (Td2(t) + Ts0(t) + Tp1(t)*sqr(2))/SQR(4) 
          
      CASE 12                            ! ● dsp3-a1　三角両錘　z軸方向  
         LET WFsh2 = (Td0(t) + Tp0(t))/SQR(2)
          
      CASE 13                            ! ● dsp3-a2　三角両錘　x軸方向
         LET WFsh2 = (Ts0(t) + Tp1(t)*sqr(2))/SQR(3) 
          
      CASE 14                            ! ● d2sp3　正八面体　z軸方向 
         LET WFsh2 = Td0(t)/SQR(3) + Ts0(t)/SQR(6) + Tp0(t)/SQR(2)  
         !   WFsh2 = Td2(t)/SQR(3) + Ts0(t)/SQR(6) + Tp1(t)/SQR(2)  x軸方向
          
      case 15                            ! ● dsp3-b1　正方錘　z軸方向
         LET WFsh2 = Td0(t)*SQR(3/10) + Ts0(t)/SQR(5) + Tp0(t)/SQR(2)  
          
      CASE 16                            ! ● dsp3-b2　正方錘　x軸方向 
         LET WFsh2 = Td2(t)*SQR(7/40) + Ts0(t)/SQR(5) + Tp1(t)/SQR(2) - Tp0(t)/SQR(8) 
          
      CASE 17                            ! ● pz
         LET WFsh2 = Tp0(t)
          
      CASE 18                            ! ● px
         LET WFsh2 = Tp1(t)
          
          
      CASE ELSE
      END SELECT
       
   CASE ELSE    ! GUI からの呼び出しではなく、DEF文で定義した関数を使用する場合
      LET WFsh2 = WFsh(t)
   END select
END FUNCTION



FUNCTION WF2(r,t)
! when guiflag = 1

   SELECT CASE guiflag
   CASE 1   ! GUIからの関数呼び出しである場合
      SELECT CASE QNN      ! 主量子数と方位量子数より、動径波動関数を指定する
      CASE 1
         LET WF2r = R1s(r)
      CASE 2
         SELECT CASE QNM
         CASE 1
            LET WF2r = R2s(r)
         CASE 2
            LET WF2r = R2p(r)
         END SELECT
      CASE 3
         SELECT CASE QNM
         CASE 1
            LET WF2r = R3s(r)
         CASE 2
            LET WF2r = R3p(r)
         CASE 3
            LET WF2r = R3d(r)
         END SELECT
      CASE 4
         SELECT CASE QNM
         CASE 1
            LET WF2r = R4s(r)
         CASE 2
            LET WF2r = R4p(r)
         CASE 3
            LET WF2r = R4d(r)
         CASE 4
            LET WF2r = R4f(r)
         END SELECT
      CASE ELSE
      END SELECT
       
      SELECT CASE QNM      ! 方位量子数と磁気量子数より、球面調和関数を指定する
      CASE 1 ! s軌道
         LET WF2 = WF2r * Ts0(t)
      CASE 2 ! p軌道
         SELECT CASE QNL
         CASE 0
            LET WF2 = WF2r * Tp0(t)
         CASE 1
            LET WF2 = WF2r * Tp1(t)
         END SELECT
      CASE 3 ! d軌道
         SELECT CASE QNL
         CASE 0
            LET WF2 = WF2r * Td0(t)
         CASE 1
            LET WF2 = WF2r * Td1(t)
         CASE 2
            LET WF2 = WF2r * Td2(t)
         END SELECT
      CASE 4 ! f軌道
         SELECT CASE QNL
         CASE 0
            LET WF2 = WF2r * Tf0(t)
         CASE 1
            LET WF2 = WF2r * Tf1(t)
         CASE 2
            LET WF2 = WF2r * Tf2(t)
         CASE 3
            LET WF2 = WF2r * Tf3(t)
         END SELECT
      CASE ELSE
      END SELECT
   CASE 2      ! GUI2からの関数呼び出しである場合 ●●●●●
      SELECT CASE FnFlag
       
      CASE 1                             ! ● 混成 p
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r = R2p(r)
         CASE 3
            LET WF2r = R3p(r)  
         CASE 4
            LET WF2r = R4p(r)
         CASE ELSE
         END SELECT
         LET WF2 = sqr(FactorA/100)*WF2r*Tp0(t) + sqr((100-FactorA)/100)*WF2r*Tp1(t)
          
      CASE 2                             ! ● sp
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = (WF2r2*Ts0(t) - WF2r*Tp0(t))/SQR(2)
          
      CASE 3                             ! ● sp2-1
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r2*Ts0(t)/SQR(3)  - WF2r*Tp0(t)*SQR(2/3)
          
      CASE 4                             ! ● sp2-2
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r2*Ts0(t)/SQR(3)  + WF2r*Tp0(t)/SQR(6) - WF2r*Tp1(t)/SQR(2)           
          
      CASE 5                             ! ● sp2-3
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r2*Ts0(t)/SQR(3)  + WF2r*Tp0(t)/SQR(6) + WF2r*Tp1(t)/SQR(2)
          
      CASE 6                             ! ● sp3-1
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r2*Ts0(t)/SQR(4)  - WF2r*Tp0(t)*SQR(3/4)       
          
          
      CASE 7                             ! ● sp3-2
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = (WF2r2*Ts0(t)*SQR(3) + WF2r*Tp0(t) - WF2r*Tp1(t)*sqr(8))/SQR(12)         
          
      CASE 8                             ! ● sp2-3
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = (WF2r2*Ts0(t)*SQR(3) + WF2r*Tp0(t) + WF2r*Tp1(t)*sqr(2))/SQR(12) 
          
      CASE 9                             ! ● sp5-1
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = (WF2r2*Ts0(t) - WF2r*Tp0(t)*SQR(3) - WF2r*Tp1(t)*sqr(2))/SQR(6) 
          
      CASE 10                            ! ● sp5-2
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r  = R2p(r)
            Let WF2r2 = R2s(r)
         CASE 3
            LET WF2r  = R3p(r)
            Let WF2r2 = R3s(r)
         CASE 4
            LET WF2r  = R4p(r)
            Let WF2r2 = R4s(r)
         CASE ELSE
         END SELECT
         LET WF2 = (WF2r2*Ts0(t) + WF2r*Tp0(t)*SQR(3) - WF2r*Tp1(t)*sqr(2))/SQR(6) 
          
      CASE 11                            ! ● dsp2　 3d2 4s 4p1  正方平面型
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r)
         LET WF2 = (WF2r3*Td2(t) + WF2r2*Ts0(t) - WF2r*Tp1(t)*sqr(2))/SQR(4) 
          
      CASE 12                            ! ● dsp3-a1　三角両錘　z軸方向
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r)
         LET WF2 = (WF2r3*Td0(t) - WF2r*Tp0(t))/SQR(2)
          
      CASE 13                            ! ● dsp3-a2　三角両錘　x軸方向
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r)
         LET WF2 = (WF2r2*Ts0(t) - WF2r*Tp1(t)*sqr(2))/SQR(3) 
          
      CASE 14                            ! ● d2sp3　正八面体　z軸方向
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r) 
         LET  WF2 = WF2r3*Td0(t)/SQR(3) + WF2r2*Ts0(t)/SQR(6) - WF2r*Tp0(t)/SQR(2)  
         !LET WF2 = WF2r3*Td2(t)/SQR(3) + WF2r2*Ts0(t)/SQR(6) - WF2r*Tp1(t)/SQR(2)  x軸方向
          
      case 15                            ! ● dsp3-b1　正方錘　z軸方向
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r)
         LET WF2 = WF2r3*Td0(t)*SQR(3/10) + WF2r2*Ts0(t)/SQR(5) - WF2r*Tp0(t)/SQR(2)  
          
      CASE 16                            ! ● dsp3-b2　正方錘　x軸方向
         LET WF2r  = R4p(r)
         Let WF2r2 = R4s(r)
         LET WF2r3 = R3d(r)
         LET WF2 = WF2r3*Td2(t)*SQR(7/40) + WF2r2*Ts0(t)/SQR(5) - WF2r*Tp1(t)/SQR(2) + WF2r*Tp0(t)/SQR(8)          
          
      CASE 17                            ! ● pz
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r = R2p(r)
         CASE 3
            LET WF2r = R3p(r)  
         CASE 4
            LET WF2r = R4p(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r*Tp0(t)
          
      CASE 18                            ! ● px
         SELECT CASE QNNhyp
         CASE 2
            LET WF2r = R2p(r)
         CASE 3
            LET WF2r = R3p(r)  
         CASE 4
            LET WF2r = R4p(r)
         CASE ELSE
         END SELECT
         LET WF2 = WF2r*Tp1(t)
          
          
      CASE ELSE                          ! ●
      END select
   CASE ELSE    ! when GUIflag = 0
      LET WF2 = WF(r,t)   ! DEF指定の関数を引用する
   END select    
end function

! ============================================================================
! GUI 表示させる（サブルーチン内でループ完結させている）。地の部分での定義は無視。
SUB GUI
   IF CautionOfGUIFlag = 0 THEN
      CALL CautionOfGUI
      LET CautionOfGUIFlag = 1
   END IF
    
   LET GUIFlag = 1   ! 関数定義でも使用　0はGUI不使用。1, 2 はそれぞれパターンA, B 
    
    
   ! 初期パラメータ等
   LET WIni = 40
   IF W = 0 THEN LET W = WIni   
   IF AOA = 0 AND AOB = 0 AND AOC = 0 AND AOD = 0 THEN LET AOD = 1
   IF SPA = 0 AND SPB = 0 THEN LET SPC = 1
   IF QNN = 0 THEN LET QNN = 4
   IF QNM = 0 THEN LET QNM = 4
   LET H03 = 1                   ! sp2-1
   LET H04 = 1                   ! sp2-2
   LET H05 = 1                   ! sp2-3
   LET FactorA = 100             ! 混成 p 用
   LET QNNHyp = 4                ! 混成用、主量子数
    
   ! 初期画面表示用セッティング
   ASK bitmap SIZE a,b           ! ASK PIXEL SIZE (0 , 0 ; 1 , 1) a,b
   LET corf = (a+1)/(b+1)        ! LET corf = a/b   , 画面表示のアスペクト比
   SET WINDOW -40*corf+40,40*corf+40,0,80
   SET TEXT FONT "ＭＳ ゴシック"
   SET TEXT HEIGHT 4
   CALL GUITextA         ! 初期画面の印字（パターンＡ）： GUIFlag = 1
   CALL GUIButtonA       ! パラメータ位置を反映
    
   ! ============================================================================
   DO
   ! GUI画面と描画画面の切り替えごとに必要
      SET WINDOW -40*corf+40,40*corf+40,0,80
      SET TEXT FONT "ＭＳ ゴシック"
      SET TEXT HEIGHT 4 
      IF GUIFlag = 1 THEN
         CALL GUITextA         ! 初期画面の印字（パターンＡ）： GUIFlag = 1
         CALL GUIButtonA       ! パラメータ位置を反映    
      ELSE
         CALL GUITextB         ! 初期画面の印字（パターンＢ）： GUIFlag = 2
         CALL GUIButtonB       ! パラメータ位置を反映 
      END IF
       
      DO   
      ! マウスを押して離すのを待ち、その位置に応じてパラメータを変える。（パターンＡ）： GUIFlag = 1
         CALL MouseCheck    ! マウスが押されて離されるまで、プログラムは待機
         IF GUIFlag = 1 THEN
            CALL GUIParamA     ! パターンAのときのパラメータの変更、マウス位置情報に応じて
            CALL GUIButtonA    ! パターンAのときのGUIボタンの再表示、パラメータに応じて
         ELSE
            CALL GUIParamB     ! パターンAのときのパラメータの変更、マウス位置情報に応じて
            CALL GUIButtonB    ! パターンAのときのGUIボタンの再表示、パラメータに応じて
         END IF        
          
         IF 75 < MouseY AND MouseY < 80 AND 45 < MouseX AND MouseX < 75  THEN ! 【画面切替】
            IF  GUIFlag = 1 THEN
               LET  GUIFlag = 2
               CALL GUITextB      ! パターンBのときの初期画面の再表示
               CALL GUIButtonB    ! パターンBのときのGUIボタンの再表示、パラメータに応じて
            ELSE
               LET  GUIFlag = 1
               CALL GUITextA      ! パターンAのときの初期画面の再表示
               CALL GUIButtonA    ! パターンAのときのGUIボタンの再表示、パラメータに応じて
            END IF
         END IF
          
         IF 0 < MouseY AND MouseY < 5 AND 45 < MouseX AND MouseX < 75 THEN !【 描画開始 】
            IF (GUIFlag = 1 AND FGUI = 1 AND SPA+SPB+SPC > 0) OR (GUIFlag = 1 AND FGUI = 0 AND AOA+AOB+AOC+AOD > 0) THEN EXIT DO
            IF  GUIFlag = 2 AND H01+H02+H03+H04+H05+H06+H07+H08+H09+H10+H11+H12+H13+H14+H15+H16+H17+H18 > 0 THEN EXIT DO
         END IF
      LOOP 
       
      ! ======================================================================================== 
      ! 【 描画開始 】選択したサブルーチンの実行
      IF GUIFlag = 1 THEN                                    ! 基底関数（混成前）, GUIFlag = 1
         LET FnFlag = 0 ! GUIFlag = 2 の時、軌道を選択するパラメータ
         LET WFSM = 0   ! 軌道種類替えるごとに初期化が必要
         CLEAR                 ! 画面初期化
         SET TEXT HEIGHT w/20  ! 画面初期化
         ! 関数を決定して、サブルーチンを実行
         IF  FGUI = 1 THEN   ! 動径波動関数の表示  WFSh2 
            PRINT 
            PRINT "球面調和関数を表示します。→ "
            PRINT  T$&"-"&STR$(QNL)
            IF SPA = 1 THEN CALL SHSqShape
            IF SPB = 1 THEN CALL SHShape
            IF SPC = 1 THEN CALL SHonSphere
         ELSE                ! 原子軌道の表示  WF
            PRINT 
            PRINT "軌道関数（動径波動関数×球面調和関数）を表示します。→ "
            PRINT  STR$(QNN)&T$&"-"&STR$(QNL)
            IF AOA = 1 THEN CALL latice
            IF AOB = 1 THEN CALL MonteCarlo
            IF AOC = 1 THEN CALL contourBasedOnInteg90
            IF AOD = 1 THEN CALL contourBasedOnInteg
         END IF    
      ELSE                                                     ! 混成関数, GUIFlag = 2
         SELECT CASE HYP
         CASE 0
            PRINT 
            PRINT "球面調和関数（二乗）を表示します。→ "
         CASE 1
            PRINT 
            PRINT "軌道関数の概形（等高線）を表示します。→ "
         CASE 2
            PRINT 
            PRINT "軌道関数の概形（点描）を表示します。→ "
         CASE ELSE
         END select
         CLEAR                  ! 画面初期化
         SET TEXT HEIGHT w/20   ! 画面初期化
         LET WFSM = 0           ! 本来、軌道種類替えるごとに初期化が必要だが z 軸に回転対称ではない軌道では、あえて、軌道毎には初期化しない。
         LET FnFlag = 0         ! ●●●●● 以下、選択された関数の表示を実行、関数定義（２か所）に書き加えること。
          
         IF H02 = 1 THEN      ! ● sp                                 !
            LET WFSM = 0      ! 軌道種類替えるごとに初期化が必要
            LET FnFlag = 2    ! 2つの関数に引き渡す
            PRINT "・ sp : 1/√2 s + 1/√2 pz"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H03 = 1 THEN      ! ● sp2-1                                 !
            LET WFSM = 0      ! z 軸回転対称なので、初期化。
            LET FnFlag = 3    ! 2つの関数に引き渡す
            PRINT "・ sp2-1 : 1/√3 s + √(2/3) pz"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H04 = 1 THEN      ! ● sp2-2                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 4    ! 2つの関数に引き渡す
            PRINT "・ sp2-2 : 1/√3 s - 1/√6 pz + 1/√2 px"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H05 = 1 THEN      ! ● sp2-3                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 5    ! 2つの関数に引き渡す
            PRINT "・ sp2-3 : 1/√3 s - 1/√6 pz - 1/√2 px"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H06 = 1 THEN      ! ● sp3-1                                 !
            LET WFSM = 0      ! z 軸回転対称なので、初期化。
            LET FnFlag = 6    ! 2つの関数に引き渡す
            PRINT "・ sp3-1 : 1/2 s + √3/2 pz"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H07 = 1 THEN      ! ● sp3-2                                !
         !  LET WFSM = 0      ! 初期かしない
            LET FnFlag = 7    ! 2つの関数に引き渡す
            PRINT "・ sp3-2 : 1/2 s - 1/√6 pz + 1/√3 px"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H08 = 1 THEN      ! ● sp3-3                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 8    ! 2つの関数に引き渡す
            PRINT "・ sp3-3 : 1/2 s - 1/√6 pz - 1/√6 px (± 1/√2 py)"
            PRINT "           ただし、py 成分は表示に反映されていません"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H09 = 1 THEN      ! ● sp5-1                                 !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 9    ! 2つの関数に引き渡す
            PRINT "・ sp5-1 : 1/√6 s + 1/√3 px  + 1/√2 pz"
            print "           ただし、sp2 が 1/√3 s - 1/√6 px  ± 1/√2 py のとき"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H10 = 1 THEN      ! ● sp5-2                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 10   ! 2つの関数に引き渡す
            PRINT "・ sp5-2 : 1/√6 s + 1/√3 px  - 1/√2 pz"
            print "           ただし、sp2 が 1/√3 s - 1/√6 px  ± 1/√2 py のとき"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H11 = 1 THEN      ! ● dsp2　 3d2 4s 4p1  正方平面型                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 11   ! 2つの関数に引き渡す
            PRINT "・ dsp2 : 1/2 3d(x2-y2)  + 1/2 4s  ± 1/√2 4px （と仮定）"
            print "           x,y 平面内、正方平面型配位、x軸方向"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H12 = 1 THEN      ! ● dsp3-a1　三角両錘　z軸方向                                !
            LET WFSM = 0      ! z軸回転対称なので、初期化する。
            LET FnFlag = 12   ! 2つの関数に引き渡す
            PRINT "・ dsp3-a1 : 1/√2 3d(z2)  ± 1/√2 4pz （と仮定）"
            print "           x,y 平面内正三角形、三角両錘型配位、z軸方向"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H13 = 1 THEN      ! ● dsp3-a2　三角両錘　x軸方向                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 13   ! 2つの関数に引き渡す
            PRINT "・ dsp3-a2 : 1/√3 4s  ± 2/√3 4px （と仮定）"
            print "           x,y 平面内正三角形、三角両錘型配位、x軸方向"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H14 = 1 THEN      ! ● d2sp3　正八面体　z軸方向                                !
            LET WFSM = 0      ! z軸回転対称なので、初期化する。
            LET FnFlag = 14   ! 2つの関数に引き渡す
            PRINT "・ d2sp3 : 1/√3 3d(z2) + 1/√6 4s ± 1/√2 4pz （と仮定）" 
            print "           正八面体　z軸方向"
            PRINT "  x 方向は : 1/√3 3d{√(x2-y2)2 + (xy)2} + 1/√6 4s ± 1/√2 4px （と仮定）"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H15 = 1 THEN      ! ● dsp3-b1　正方錘　z軸方向                                !
            LET WFSM = 0      ! z軸回転対称なので、初期化する。
            LET FnFlag = 15   ! 2つの関数に引き渡す
            PRINT "・ dsp3-b1 : √(3/10) 3d(z2) + 1/√5 4s + 1/√2 4pz （と仮定）"
            print "           正方錘型配位、z軸方向"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H16 = 1 THEN      ! ● dsp3-b2　正方錘　x軸方向                                !
         !  LET WFSM = 0      ! 初期化しない
            LET FnFlag = 16   ! 2つの関数に引き渡す
            PRINT "・ dsp3-b2 : √(7/40) 3d{√(x2-y2)2 + (xy)2} + 1/√5 4s + 1/√2 4px - 1/√8 Pz（と仮定）" 
            print "           正方錘型配位、x軸方向"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H17 = 1 THEN      ! ● pz
            LET WFSM = 0      ! 初期化する                                 !
            LET FnFlag = 17   ! 2つの関数に引き渡す
            PRINT "・ p 軌道：pz"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H18 = 1 THEN      ! ● px    
         !  LET WFSM = 0      ! 初期化しない                             !
            LET FnFlag = 18   ! 2つの関数に引き渡す
            PRINT "・ p 軌道：px"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
         IF H01 = 1 THEN      ! ● 混成 p                                 !
         !  LET WFSM = 0      ! 初期化しない                             !
            LET FnFlag = 1    ! 2つの関数に引き渡す
            PRINT "・ 混成 p 軌道：pz"; FactorA; "%, px"; 100-FactorA ; "% の寄与"
            IF Hyp = 0 THEN   CALL SHSqShape
            IF Hyp = 1 THEN   CALL contourBasedOnWFSM
            IF Hyp = 2 THEN   CALL MonteCarlo
         END IF
          
      END IF
       
      PRINT 
      print "指定した表示は終了しました。"
      pause "「OK 続行」で GUI に戻る（図は消えます）。"& CHR$(13) &"「×で閉じる」と、デバッグモード、図をコピーできます。"& CHR$(13) &"その後「続行」でGUIに戻ります。"
      ! STOP ! ← 最後にトル
   loop
END SUB

SUB GUITextA
! 初期画面の印字（パターンＡ）： GUIFlag = 1
   CLEAR
   PLOT TEXT ,AT 0, 75 : " 基底関数"
   PLOT TEXT ,AT 45, 75 : "【 画面切替 】"
   PLOT TEXT ,AT 0, 70 : "○ 球面調和関数の表示"        ! FGUI = 1
   PLOT TEXT ,AT 0, 65 : "　　□ 絶対値"                   ! SPA = 1/0
   PLOT TEXT ,AT 0, 60 : "　　□ 二乗"                     ! SPB = 1/0
   PLOT TEXT ,AT 0, 55 : "　　□ 円周上の小円"             ! SPC = 1/0
   PLOT TEXT ,AT 0, 50 : "○ 原子軌道の表示         "   ! FGUI = 0
   PLOT TEXT ,AT 0, 45 : "　　□ 格子点"                   ! AOA = 1/0
   PLOT TEXT ,AT 0, 40 : "　　□ 点描法"                   ! AOB = 1/0
   PLOT TEXT ,AT 0, 35 : "　　□ 等高線（90）"             ! AOC = 1/0
   PLOT TEXT ,AT 0, 30 : "　　□ 等高線（10〜90）"         ! AOD = 1/0
   PLOT TEXT ,AT 0, 25 : "軌道の選択（上から順に選択↓）"
   PLOT TEXT ,AT 0, 20 : "→ 主量子数   [   ] ："          ! QNN
   PLOT TEXT ,AT 55, 20 : "1,2,3,4"
   PLOT TEXT ,AT 0, 15 : "→ 方位量子数 [   ] ："          ! QNM 1,2,3,4
   PLOT TEXT ,AT 55, 15 : "s,p,d,f"
   PLOT TEXT ,AT 0, 10 : "→ 磁気量子数 [   ] ："          ! QNL
   PLOT TEXT ,AT 55, 10 : "0,1,2,3"
   PLOT TEXT ,AT 0,  5 : "→ 画面サイズ [   ] ："          ! w  10 〜 60 step 5
   PLOT TEXT ,AT 55, 5 : "≪ ･ ≫"
   PLOT TEXT ,AT 0,  0 : "　　　　　　　　【 描画開始 】"
END SUB

SUB GUIButtonA   ! パターンAのときのGUIボタンの再表示、パラメータに応じて
   SET TEXT COLOR 0
   PLOT TEXT ,AT 0, 70 : "●"
   PLOT TEXT ,AT 0, 65 : "　　■"
   PLOT TEXT ,AT 0, 60 : "　　■" 
   PLOT TEXT ,AT 0, 55 : "　　■"
   PLOT TEXT ,AT 0, 50 : "●"
   PLOT TEXT ,AT 0, 45 : "　　■"
   PLOT TEXT ,AT 0, 40 : "　　■"
   PLOT TEXT ,AT 0, 35 : "　　■"
   PLOT TEXT ,AT 0, 30 : "　　■"
   SET TEXT HEIGHT 6
   PLOT TEXT ,AT 38, 19 : "■"
   PLOT TEXT ,AT 38, 14 : "■"
   PLOT TEXT ,AT 38,  9 : "■"
   PLOT TEXT ,AT 38,  4 : "■"
   SET TEXT HEIGHT 4
   SET TEXT COLOR 1
   IF FGUI = 1 THEN
      PLOT TEXT ,AT 0, 70 : "●"
      PLOT TEXT ,AT 0, 50 : "○"
   ELSE
      PLOT TEXT ,AT 0, 70 : "○"
      PLOT TEXT ,AT 0, 50 : "●"
   END IF
   IF FGUI = 1 AND SPA = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 65 : "　　" & SQ$
   IF FGUI = 1 AND SPB = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 60 : "　　" & SQ$
   IF FGUI = 1 AND SPC = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 55 : "　　" & SQ$
   IF FGUI = 0 AND AOA = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 45 : "　　" & SQ$
   IF FGUI = 0 AND AOB = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 40 : "　　" & SQ$
   IF FGUI = 0 AND AOC = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 35 : "　　" & SQ$
   IF FGUI = 0 AND AOD = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 0, 30 : "　　" & SQ$
    
   PLOT TEXT ,AT 39, 20 : " " & STR$(QNN)                                            ! 主量子数 
   SELECT CASE QNM
   CASE 1
      LET T$ ="s"
   CASE 2
      LET T$ ="p"
   CASE 3
      LET T$ ="d"
   CASE 4
      LET T$ ="f"
   END SELECT
   PLOT TEXT ,AT 39, 15 : " " & T$
   PLOT TEXT ,AT 39, 10 : " " & STR$(QNL)
   IF FGUI = 0 THEN PLOT TEXT ,AT 39,  5 : REPEAT$(" ",2-LEN(STR$(w))) & STR$(w)    ! 画面サイズ
    
END SUB

SUB GUIParamA  ! パターンAのときのパラメータの変更、マウス位置情報に応じて
   IF 70 < MouseY AND MouseY < 75 AND 0 < MouseX AND MouseX < 5 THEN LET FGUI = 1
   IF 50 < MouseY AND MouseY < 55 AND 0 < MouseX AND MouseX < 5 THEN LET FGUI = 0
   IF FGUI = 1 AND 66 < MouseY AND MouseY < 70 AND 10 < MouseX AND MouseX < 15 THEN 
      IF SPA = 0 THEN LET SPA = 1 ELSE LET SPA = 0
   END IF
   IF FGUI = 1 AND 61 < MouseY AND MouseY < 65 AND 10 < MouseX AND MouseX < 15 THEN
      IF SPB = 0 THEN LET SPB = 1 ELSE LET SPB = 0
   END IF          
   IF FGUI = 1 AND 56 < MouseY AND MouseY < 60 AND 10 < MouseX AND MouseX < 15 THEN
      IF SPC = 0 THEN LET SPC = 1 ELSE LET SPC = 0
   END IF  
   IF FGUI = 0 AND 46 < MouseY AND MouseY < 50 AND 10 < MouseX AND MouseX < 15 THEN
      IF AOA = 0 THEN
         LET AOA = 1
         LET AOB = 0
      ELSE
         LET AOA = 0
      END IF 
   END IF
   IF FGUI = 0 AND 41 < MouseY AND MouseY < 45 AND 10 < MouseX AND MouseX < 15 THEN
      IF AOB = 0 THEN
         LET AOB = 1
         LET AOA = 0
      ELSE
         LET AOB = 0
      END IF 
   END IF             
   IF FGUI = 0 AND 36 < MouseY AND MouseY < 40 AND 10 < MouseX AND MouseX < 15 THEN
      IF AOC = 0 THEN
         LET AOC = 1
         LET AOD = 0
      ELSE
         LET AOC = 0
      END IF 
   END IF
   IF FGUI = 0 AND 31 < MouseY AND MouseY < 35 AND 10 < MouseX AND MouseX < 15 THEN
      IF AOD = 0 THEN
         LET AOD = 1
         LET AOC = 0
      ELSE
         LET AOD = 0
      END IF 
   END IF     
   IF 21 < MouseY AND MouseY < 25 THEN ! 主量子数
      SELECT CASE MouseX
      CASE 54 TO 59
         LET QNN = 1
         LET QNM = 1
         LET QNL = 0
      CASE 59 TO 64
         LET QNN = 2
         IF QNM = 3 OR QNM = 4 THEN 
            LET QNM = 1
            LET QNL = 0
         END IF
      CASE 64 TO 69
         LET QNN = 3
         IF QNM = 4 THEN
            LET QNM = 1
            LET QNL = 0
         END IF
      CASE 69 TO 74
         LET QNN = 4
      CASE ELSE
      END SELECT
   END IF
   IF 16 < MouseY AND MouseY < 20 THEN ! 方位量子数
      SELECT CASE MouseX
      CASE 54 TO 59
         LET QNM = 1
         LET QNL = 0
      CASE 59 TO 64
         IF QNN => 2 THEN
            LET QNM = 2
            IF QNL > 2 THEN LET QNL = 1
         END IF
      CASE 64 TO 69
         IF QNN => 3 THEN
            LET QNM = 3
            IF QNL = 4 THEN LET QNL = 1
         END IF
      CASE 69 TO 74
         IF QNN => 4 THEN LET QNM = 4
      CASE ELSE
      END SELECT
   END IF
   IF 11 < MouseY AND MouseY < 15 THEN ! 磁気量子数
      SELECT CASE MouseX 
      CASE 54 TO 59
         LET QNL = 0
      CASE 59 TO 64
         IF QNM => 2 THEN LET QNL = 1
      CASE 64 TO 69
         IF QNM => 3 THEN LET QNL = 2
      CASE 69 TO 74
         IF QNM => 4 THEN LET QNL = 3
      CASE ELSE
      END SELECT
   END IF          
   IF 6 < MouseY AND MouseY < 10 AND FGUI = 0 THEN ! 画面サイズ w
      SELECT CASE MouseX
      CASE 54 TO 62
         LET W = W - 5
         IF W < 5 THEN LET W = 5
      CASE 62 TO 66
         LET W = 40
      CASE 66 TO 74
         LET W = W + 5
         IF W > 60 THEN LET W = 60
      CASE ELSE
      END SELECT
   END IF              
END SUB

SUB GUITextB
! 初期画面の印字（パターンＢ）： GUIFlag = 2
   CLEAR
   PLOT TEXT ,AT 0, 75 : " 混成関数"
   PLOT TEXT ,AT 45, 75 : "【 画面切替 】"
    
   PLOT TEXT ,AT 0, 70 : "○ 球面調和関数の二乗"         ! HYP = 0
   PLOT TEXT ,AT 0, 65 : "○ 軌道の概形（等高線）"       ! HYP = 1  
   PLOT TEXT ,AT 0, 60 : "○ 軌道の概形（点描）"         ! HYP = 2
   PLOT TEXT ,AT 5, 55 : "□ 混成 p ： √A･pz + √B･px"  ! H01 = 1/0
   PLOT TEXT ,AT 30, 50 : "A%[    ] "                    ! FactorA %
   PLOT TEXT ,AT 55, 50 : "≪ ･ ≫"
   PLOT TEXT ,AT  5, 45 : "□ sp"                        ! H02 = 1/0
   PLOT TEXT ,AT 30, 45 : "□ pz(B=0)"                     ! H17 = 1/0
   PLOT TEXT ,AT 55, 45 : "□ px(A=0)"                     ! H18 = 1/0
   PLOT TEXT ,AT  5, 40 : "□ sp2-1"                     ! H03 = 1/0
   PLOT TEXT ,AT 30, 40 : "□ sp2-2"                     ! H04 = 1/0
   PLOT TEXT ,AT 55, 40 : "□ sp2-3"                     ! H05 = 1/0
   PLOT TEXT ,AT  5, 35 : "□ sp3-1"                     ! H06 = 1/0
   PLOT TEXT ,AT 30, 35 : "□ sp3-2"                     ! H07 = 1/0
   PLOT TEXT ,AT 55, 35 : "□ sp3-3"                     ! H08 = 1/0
   PLOT TEXT ,AT  5, 30 : "□ sp5-1"                     ! H09 = 1/0
   PLOT TEXT ,AT 30, 30 : "□ sp5-2"                     ! H10 = 1/0
   !           
   SET TEXT COLOR 15
   PLOT TEXT ,AT  5, 23 : "□ dsp2"                      ! H11 = 1/0 
   PLOT TEXT ,AT 30, 23 : "□ dsp3-a1"                   ! H12 = 1/0
   PLOT TEXT ,AT 55, 23 : "□ dsp3-a2"                   ! H13 = 1/0   
   PLOT TEXT ,AT  5, 18 : "□ d2sp3"                     ! H14 = 1/0 
   PLOT TEXT ,AT 30, 18 : "□ dsp3-b1"                   ! H15 = 1/0    
   PLOT TEXT ,AT 55, 18 : "□ dsp3-b2"                   ! H16 = 1/0      
    
   SET TEXT COLOR 1
   PLOT TEXT ,AT 0, 10 : "→ 主量子数   [   ] ："        ! QNNHyp
   PLOT TEXT ,AT 60, 10 : "2,3,4"
   PLOT TEXT ,AT 0,  5 : "→ 画面サイズ [   ] ："        ! w  10 〜 60 step 5
   PLOT TEXT ,AT 55, 5 : "≪ ･ ≫"
   PLOT TEXT ,AT 0,  0 : "　　　　　　　　【 描画開始 】"
END SUB

SUB GUIButtonB
   SET TEXT COLOR 0
   PLOT TEXT ,AT  0, 70 : "●"   
   PLOT TEXT ,AT  0, 65 : "●"     
   PLOT TEXT ,AT  0, 60 : "●" 
   PLOT TEXT ,AT  5, 55 : "■"
   PLOT TEXT ,AT  5, 45 : "■"
   PLOT TEXT ,AT 30, 45 : "■"
   PLOT TEXT ,AT 55, 45 : "■"
   PLOT TEXT ,AT  5, 40 : "■"
   PLOT TEXT ,AT 30, 40 : "■"
   PLOT TEXT ,AT 55, 40 : "■"
   PLOT TEXT ,AT  5, 35 : "■"
   PLOT TEXT ,AT 30, 35 : "■"
   PLOT TEXT ,AT 55, 35 : "■"
   PLOT TEXT ,AT  5, 30 : "■"
   PLOT TEXT ,AT 30, 30 : "■"
   PLOT TEXT ,AT  5, 23 : "■"
   PLOT TEXT ,AT 30, 23 : "■"
   PLOT TEXT ,AT 55, 23 : "■"
   PLOT TEXT ,AT  5, 18 : "■"
   PLOT TEXT ,AT 30, 18 : "■"
   PLOT TEXT ,AT 55, 18 : "■"
   SET TEXT HEIGHT 6  
   PLOT TEXT ,AT 38, 49 : "■"         ! 係数A 
   PLOT TEXT ,AT 40, 49 : "■"         ! 係数A 
   PLOT TEXT ,AT 38,  9 : "■"         ! 主量子数  
   PLOT TEXT ,AT 38,  4 : "■"
    
   SET TEXT COLOR 1  
   SET TEXT HEIGHT 4
    
   SELECT CASE HYP
   CASE 0
      PLOT TEXT ,AT 0, 70 : "●"       ! HYP = 0
      PLOT TEXT ,AT 0, 65 : "○"       ! HYP = 1  
      PLOT TEXT ,AT 0, 60 : "○"       ! HYP = 2
   CASE 1
      PLOT TEXT ,AT 0, 70 : "○"       ! HYP = 0
      PLOT TEXT ,AT 0, 65 : "●"       ! HYP = 1  
      PLOT TEXT ,AT 0, 60 : "○"       ! HYP = 2
   CASE 2
      PLOT TEXT ,AT 0, 70 : "○"       ! HYP = 0
      PLOT TEXT ,AT 0, 65 : "○"       ! HYP = 1  
      PLOT TEXT ,AT 0, 60 : "●"       ! HYP = 2
   CASE ELSE
   END select
    
   IF H01 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 55 : SQ$
   IF H02 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 45 : SQ$
   IF H03 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 40 : SQ$
   IF H04 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 40 : SQ$
   IF H05 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 55, 40 : SQ$
   IF H06 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 35 : SQ$
   IF H07 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 35 : SQ$
   IF H08 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 55, 35 : SQ$
   IF H09 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 30 : SQ$
   IF H10 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 30 : SQ$
   IF H11 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 23 : SQ$
   IF H12 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 23 : SQ$
   IF H13 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 55, 23 : SQ$
   IF H14 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT  5, 18 : SQ$
   IF H15 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 18 : SQ$
   IF H16 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 55, 18 : SQ$
   IF H17 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 30, 45 : SQ$
   IF H18 = 1 THEN LET SQ$ = "■" ELSE LET SQ$ = "□"
   PLOT TEXT ,AT 55, 45 : SQ$ 
   IF HYP <> 0 THEN PLOT TEXT ,AT 39, 10 : " " & STR$(QNNHyp)                                         ! 主量子数  
   IF HYP <> 0 THEN PLOT TEXT ,AT 39,  5 : REPEAT$(" ",2-LEN(STR$(w))) & STR$(w)                      ! 画面サイズ
   IF H01 = 1  THEN PLOT TEXT ,AT 39, 50 : REPEAT$(" ",3-LEN(STR$(FactorA))) & STR$(FactorA)      ! 係数A 
    
END SUB


SUB GUIParamB

   IF 70 < MouseY AND MouseY < 75 AND 0 < MouseX AND MouseX < 5 THEN LET HYP = 0
   IF 65 < MouseY AND MouseY < 70 AND 0 < MouseX AND MouseX < 5 THEN LET HYP = 1
   IF 60 < MouseY AND MouseY < 65 AND 0 < MouseX AND MouseX < 5 THEN LET HYP = 2
    
   IF 55 < MouseY AND MouseY < 60 AND  5 < MouseX AND MouseX < 10 THEN LET H01 = MOD(H01+1,2)
   IF 45 < MouseY AND MouseY < 50 AND  5 < MouseX AND MouseX < 10 THEN LET H02 = MOD(H02+1,2)
   IF 40 < MouseY AND MouseY < 45 AND  5 < MouseX AND MouseX < 10 THEN LET H03 = MOD(H03+1,2)
   IF 40 < MouseY AND MouseY < 45 AND 30 < MouseX AND MouseX < 35 THEN LET H04 = MOD(H04+1,2)
   IF 40 < MouseY AND MouseY < 45 AND 55 < MouseX AND MouseX < 60 THEN LET H05 = MOD(H05+1,2)
   IF 35 < MouseY AND MouseY < 40 AND  5 < MouseX AND MouseX < 10 THEN LET H06 = MOD(H06+1,2)
   IF 35 < MouseY AND MouseY < 40 AND 30 < MouseX AND MouseX < 35 THEN LET H07 = MOD(H07+1,2)
   IF 35 < MouseY AND MouseY < 40 AND 55 < MouseX AND MouseX < 60 THEN LET H08 = MOD(H08+1,2)
   IF 30 < MouseY AND MouseY < 35 AND  5 < MouseX AND MouseX < 10 THEN LET H09 = MOD(H09+1,2)
   IF 30 < MouseY AND MouseY < 35 AND 30 < MouseX AND MouseX < 35 THEN LET H10 = MOD(H10+1,2)
    
   IF 23 < MouseY AND MouseY < 28 AND  5 < MouseX AND MouseX < 10 THEN LET H11 = MOD(H11+1,2)
   IF 23 < MouseY AND MouseY < 28 AND 30 < MouseX AND MouseX < 35 THEN LET H12 = MOD(H12+1,2)
   IF 23 < MouseY AND MouseY < 28 AND 55 < MouseX AND MouseX < 60 THEN LET H13 = MOD(H13+1,2)
   IF 18 < MouseY AND MouseY < 23 AND  5 < MouseX AND MouseX < 10 THEN LET H14 = MOD(H14+1,2)
   IF 18 < MouseY AND MouseY < 23 AND 30 < MouseX AND MouseX < 35 THEN LET H15 = MOD(H15+1,2)
   IF 18 < MouseY AND MouseY < 23 AND 55 < MouseX AND MouseX < 60 THEN LET H16 = MOD(H16+1,2)
    
   IF 45 < MouseY AND MouseY < 50 AND 30 < MouseX AND MouseX < 35 THEN LET H17 = MOD(H17+1,2)
   IF 45 < MouseY AND MouseY < 50 AND 55 < MouseX AND MouseX < 60 THEN LET H18 = MOD(H18+1,2)
    
   IF 50 < MouseY AND MouseY < 55 AND H01 = 1 THEN       ! FactorA
      SELECT CASE MouseX
      CASE 54 TO 62
         LET FactorA = FactorA - 5
         IF FactorA < 0 THEN LET FactorA = 0
      CASE 62 TO 66
         LET FactorA = 50
      CASE 66 TO 74
         LET FactorA = FactorA + 5
         IF FactorA > 100 THEN LET FactorA = 100
      CASE ELSE
      END SELECT
   END IF              
    
   IF 10 < MouseY AND MouseY < 15 AND HYP <> 0 THEN    ! 主量子数
      SELECT CASE MouseX
      CASE 59 TO 64
         LET QNNHyp = 2
      CASE 64 TO 69
         LET QNNHyp = 3
      CASE 69 TO 74
         LET QNNHyp = 4
      CASE ELSE
      END SELECT
   END IF 
   ! d 軌道を含む混成を一つでも選んでいたら主量子数を 3 に強制変更     
   ! IF H11 + H12 + H13 + H14 + H15 + H16 > 0 THEN LET QNNHyp = 3
    
   IF 6 < MouseY AND MouseY < 10 AND Hyp<>0 THEN       ! 画面サイズ w
      SELECT CASE MouseX
      CASE 54 TO 62
         LET W = W - 5
         IF W < 5 THEN LET W = 5
      CASE 62 TO 66
         LET W = 40
      CASE 66 TO 74
         LET W = W + 5
         IF W > 60 THEN LET W = 60
      CASE ELSE
      END SELECT
   END IF                  
    
END SUB

SUB MouseCheck ! マウスボタンのクリックを待ち、マウス位置情報を得る
! マウスボタンが押されるまで待つ
   DO
      MOUSE POLL MouseX, MouseY, MouseLeft, MouseRight
   LOOP WHILE MouseLeft = 0
   ! マウスボタンが離れるまで待つ
   DO
      MOUSE POLL MouseX, MouseY, MouseLeft, MouseRight
   LOOP WHILE MouseLeft = 1
END SUB

SUB CautionOfGUI


   PRINT "GUI に表示されない軌道や関数の表示等については、GUI をオフにして下さい。"
   PRINT "そのためには、プログラム中「CALL GUI」の前に「!」を書き込んで下さい。"
   PRINT "グラフィックスウィンドウは、正方形かまたは縦長のサイズを選んでください。"
   PRINT
   PRINT "軌道を選択し、球面調和関数と軌道を表示できます。"
   PRINT "はじめに、○（ラジオボタン）を選んでください。"
   PRINT "続いて、表示させたい種類を、□（チェックボックス）で選んでください。"
   PRINT "[  ] の中のパラメータは、右側の領域をクリックして選択、または増減して下さい。"
   PRINT "最後に右下の【 描画開始 】ボタンを押して下さい。"
   PRINT "その都度計算していますので、等高線表示等、時間が掛かるものがあります。"
   PRINT
   PRINT "表示終了後は、ポップアップウィンドウから"
   PRINT "「OK 続行」を押すと、GUI に戻ります。"
   PRINT "表示された図を利用する場合は、"
   PRINT "「×で閉じる」と、デバッグモードに移行し、"
   PRINT "図をコピーできます。デバッグモードからは、"
   PRINT "「続行」ボタンを押せば、GUIに戻ります。"
   PRINT
   PRINT "→ 原子軌道は、z軸（表示上のy軸）に対する回転体として定義されています。"
   PRINT "たとえば、pz → p0, px + py → p1 と表示します（正確には二乗）。"
   PRINT "表示上、y 成分（手前方向）は無視し、xz 平面による断面のみです。"
   PRINT "核電荷1で計算しています。軸の数値は、ボーア半径です。"
   PRINT "→ 混成軌道は、表示した式に相当するものを表示しています。"
   PRINT "ただし、z軸に回転対称ではない軌道は、積分を正しく計算していないので、"
   PRINT "等高線表示は、初めに表示した軌道に対して、仮に得られた積分値を基に、"
   PRINT "同じ関数値になるところを結んでいますので、% 表示にはなりません。"
   PRINT "→ 特に、d 軌道を含む混成軌道は、式さえ定義できれば、表示できますが、"
   PRINT "実際に結合に寄与する軌道であることを意味するわけではありません。"
   PRINT    
END SUB
! ============================================================================

END


! 改訂履歴
! 2020.03.13 サブルーチン MonteCarlo をはじめに呼び出すと、パラメータ WFSM がゼロであるために、ゼロ除算のエラーが出る不具合を修正。
! 2020.03.18 f1の球面調和関数が間違っていました。10.5 を 11.5 と書いてしまっていましたので修正。
! 2022.01.06 基本的な原子軌道をマウスクリックで選択できるように、GUIを設定。
!            細かい設定や、混成軌道の表示などについては、GUIをoff にしてご使用下さい。
! 2022.01.18 GUI から異なる軌道を表示させる際にも WFSM 波動関数の二乗値の最大値を計算し直さない点を修正。
!            一回の起動で一種類しか描かない設定のとき、複数の等高線で計算しなおさなくてよいようにしていたが、GUI 実装時に修正忘れ。