aUCBLogo Demos and Tests / harmonic_oscillator

be harmonic_oscillator
   be psi x E
      ; this is an approximation of the real solution
      if E==0 [output 0]

;      root=sqrt -E+0i
;      tmp=(ln root)*E^2

;      root=(sqrt ((x^2)-E))
      root=(sqrt (0i+(sqr x)-E))
;setpixelXY x*100 100*real root 4
;setpixelXY x*100 100*imag root 1
;      tmp=(exp (-(abs x)*root-E*(ln -(abs x)+root+0i)-tmp)*0.5)
      tmp=(exp (-(abs x)*root/4-E*(ln 0i-(abs x)+root)/4+E*(ln 0i-E)/8)*2)
      output tmp
   end

   be expo x E
      ; this is an approximation of the real solution
      if E==0 [output 0]

      root=sqrt 0i-E
      tmp=(ln root)*sqr E

;      root=(sqrt ((x^2)-E))
      if x==0 [x=0i]
      root=(sqrt ((sqr x)-E+0i))
;setpixelXY (real x)*100 100*real root 4
;setpixelXY (real x)*100 100*imag root 1
;      tmp=(exp (-(abs x)*root-E*(ln -(abs x)+root+0i)-tmp)*0.5)
      tmp=(-(abs x)*root/4-E*(ln 0i-(abs x)+root)/4+E*(ln 0i-E)/8)*2
      output tmp
   end

   be d2phidx2 x E
      epsilon=1e-3
      y1=((sqrt psi x+epsilon E)-(sqrt psi x E))/epsilon
      y2=((sqrt psi x E)-(sqrt psi x-epsilon E))/epsilon
      output (y1-y2)/epsilon
   end

   be hermite x n
      ;this is the solution of the harmonic oscillator
      s=1
      psi0=(exp -(sqr x/s)/2)
      output case n
      [   [0  1  *psi0]
         [1  0.6*psi0*2*x]
         [2  0.25*psi0*(-4+4*x^2)]
         [3  0.2*psi0*(-12*x+8*x^3)]
         [4  0.083*psi0*(12-48*x^2+16*x^4)]
         [5   0.027*psi0*(120*x-160*x^3+32*x^5)]
         [6   0.0083*psi0*(-120+720*x^2-480*x^4+64*x^6)]
         [else 0]
      ]
   end

;   be E n
;      output 2*(n+1/2)
;   end

   be drawPsi E
      PenUp
      setXY -400 200*real psi -4 E
      PenDown
      psi0=Norm psi 0 E
      for [x -4 4 0.1]
      [   setXY x*100 100*(real psi x E)/psi0
      ]
   end

   be findRoot f x E
      y=1.5
      xo=0
      yo=10000000
      dx=1
      dy=0
      x=x+0i
      repeat 100
      [
         y=real run (List f x E)
         yx=real run (List f x+dx E)
         yy=real run (List f x+dy*1i E)
         d=sqrt (sqr dx)+sqr dy
         if d != 0
         [   dx=1e-5*(yx-y)/d
            if (Norm dx) < 0.1 [dx=0.1*dx/Norm dx]
         ]
;         if d != 0
;         [   dy=1e-5*(yy-y)/d
;            if (Norm dy) < 0.1 [dy=0.1*dy/Norm dy]
;         ]
         if y < yo
         [   yo=y
            xo=x
         ]
         x=x+dx ;+dy*1i
;setPixelXY real x imag x random IntMax
;(pr x "\    y "\    dx "\    dy)
         if (Norm x) > 1e15 [break]
      ]
(pr xo yo)
;updateGraph
      if (Norm yo) > 1e5 [yo=0]
      output yo+0i
   end

   _Z=1+0i
   C=1
   D=0
   F=-12

   be psi2 x E
   ;   output 72*(sqr _Z)*C+6*(sqr _Z)*D+F*(sqr exp _Z)-6*sqr expo x E
      output 72*(sqr _Z)*C+6*(sqr _Z)*D+F*(sqr exp _Z)-6*sqr expo x E
   end

   be drawPsi2 E
ignore [
      PenUp
      for [re -300 300 10]
      [   for [im -300 300 10]
         [   z=(re+1i*im)
            y=0.1*real psi2 z E
            setXY re im
            setFC HSB y 0.5 1
            fillRect [-5 -5][5 5]
         ]
      ]
]
      PenUp
      setXY -400 200*real psi2 -4 E
      PenDown
      tmp=(exp findRoot "psi2 -4 E)^4
      psi0=real tmp
      for [x -4 4 0.1]
      [   tmp=(exp findRoot "psi2 x E)^4
         setXY x*100 1e-10*(real tmp) ;/psi0
         if key? [break]
      ]
   end

   be drawPsi0 n
      PenUp
      setXY -400 200*real hermite -4 n
      PenDown
      for [x -4 4 0.1]
      [   setXY x*100 100*(hermite x n)
      ]
   end

   n=4
   E_=9
   setFPUStatusFlags false false false false
   hideTurtle
   setUpdateGraph false
   WindowMode
;norefresh
   E_Indicator=(FloatControl [] [E] 0 E_ 100 1 2 [] 3 [] [30 0])
   E_slider=(Slider [] [E] 1 E_*100 1000
      [   clearScreen
         clearText
         setPC 0
         E_=SliderValue/100
         FloatControlSetValue E_Indicator E_
         setPC 1
         drawPsi E_
         setPC 2
      ;   drawPsi2 E_
         setpc 4
         drawPsi0 n
         updateGraph
      ] wxSL_VERTICAL [0 0][30 400])

   clearScreen
   clearText
   setPC 1
   drawPsi E_
   setPC 2
;   drawPsi2 E_
   setpc 4
   drawPsi0 n
   updateGraph
end

be W_ x
   w=0
   repeat 100
   [   wexpw=w*exp w
      wPlusOneTimesExpW=(w+1)*exp w
      n1=2*w+2
      if n1 != 0
      [   n2=wPlusOneTimesExpW-(w+2)*(wexpw-x)/(2*w+2)
         if n2 != 0
         [   w=w-(wexpw-x)/n2
         ]
      ]
   ]
   output w
end