aUCBLogo Demos and Tests / fermisurface3


be fermisurface3
   
A1=1
   
G=2
   
a=360/200
   
Ef=0
   
kstart=-300
   
kend=300
   
dkz=10
   
dkx=10   ;when smaller then less gaps, but also slower
   
dkx_s=2   ;for the surface
   
dkz_s=2
   
dkz_s2=dkz_s*1.05   ;because the x-coordinates are different for every z
                  ;so as to compensate some ugly holes in the graphics
   
dkxfind=1
   
dkx2min=1
   
dkx2min_s=1
   
dmin=50
   
epsilon=0.001
   
branches=6
   
   
be kx ky kz
      
output -A-G*((cos kx*a)+(cos ky*a)+(cos kz*a))
   
end
   
   
be E_solved_to_ky Ef kx kz d f
;output f*sqrt 10000-sqr kx
      
tmp= -(Ef+A1)/G-(cos kx*a)-(cos kz*a)
      
ifElse tmp >= -and2 tmp <=1
      
[   output (d+f*ArcCos tmp)/a
      
][   output 99999999
      
]
   
end

   
be dE_solved_to_ky_dkx Ef kx kz d f
      
output 
      
( (E_solved_to_ky Ef kx+epsilon kz d f)
      
- (E_solved_to_ky Ef kx kz d f) )        / epsilon
   
end
   
   
be find_dE_dkx_equal_1 Ef kz d f
      
l=[]
      
for [kx kstart kend dkxfind]
      
[   tmp=abs (abs dE_solved_to_ky_dkx Ef kx kz d f)-1
         
if tmp 0.1
         
[   push "l list tmp kx
         
]
      
]
      
l=CombSort l
      
l2=[]
      
repeat branches
      
[   if l==[] [break]
         
push "l2 l.(1).(2)
         
l=butFirst l
      
]
      
output l2
   
end

   
noRefresh
   
clearText
   
setUpdateGraph false
   
clearScreen
   
WindowMode
   
hideTurtle
   
perspective
   
setPenSize [2 2]
   
disableCylinderLines
   
enableCylinderLines
   
disableRoundLineEnds
;   axes
   
   
y=0
   
ox=-100
   
oy=0
   
oz=-100
   
ok=true
   
kx=0
   
ky=0
   
kz=0
   
d=0
   
f=0
   
dkx2=0

   
be draw Ef
      
for [-360*1 360*1 360]
      
[   for [-1 1 2]
         
[   for [kz kstart kend dkz]
            
[   l=find_dE_dkx_equal_1 Ef kz d f
               
while [!= []]
               
[   ok=true
                  
kx=first l
                  
PenUp
                  
setXYZ kx kz E_solved_to_ky Ef kx kz d f
                  
dkx2=dkx
                  
while [kx <= kend]
                  
[   y=E_solved_to_ky Ef kx kz d f
                     
onePoint
                     
if == 99999999 [PenUp break]
                     
ifelse dkx2 0.001 
                     
and2 (abs (y-oy)/dkx2) > dkx2min
                     
and2 (abs y-oy) < kend
                     
[   dkx2=dkx2min*dkx2/(abs y-oy)
                     
][   dkx2=dkx
                     
]
      
ignore [
                     
(pr dkx2
                        
dkx2 0.001
                        
(abs y-oy)/dkx2 dkx2min
                        
(abs y-oy) < kend
                     
)
      
]
                     
kx=kx+dkx2
                     
oy=y
                  
]            
                  
kx=first l
                  
dkx2=dkx
                  
PenUp
                  
setXYZ kx kz E_solved_to_ky Ef kx kz d f
                  
while [kx >= kstart]
                  
[   y=E_solved_to_ky Ef kx kz d f
                     
onePoint
                     
if == 99999999 [PenUp break]
                     
ifelse dkx2 0.001 
                     
and2 (abs (y-oy)/dkx2) > dkx2min
                     
and2 (abs y-oy) < kend
                     
[   dkx2=dkx2min*dkx2/(abs y-oy)
                     
][   dkx2=dkx
                     
]
                     
kx=kx-dkx2
                     
oy=y
                  
]            
                  
l=butFirst l
               
]
            
]
         
]
      
]
;   comment [
      
for [-360*1 360*1 360]
      
[   for [-1 1 2]
         
[   PenUp
            
for [kz kstart kend dkz_s]
            
[   l=find_dE_dkx_equal_1 Ef kz d f
               
while [!= []]
               
[   ok=true
                  
kx=first l
                  
dkx2=dkx_s
                  
oy=100
                  
while [kx <= kend]
                  
[   oneQuad Ef
                     
kx=kx+dkx2
                  
]            
                  
kx=0
                  
dkx_s=-dkx_s
                  
dkx2=dkx_s
                  
oy=100
                  
while [kx >= kstart]
                  
[   oneQuad Ef
                     
kx=kx+dkx2
                  
]            
                  
dkx_s=-dkx_s
                  
l=butFirst l
               
]
            
]
         
]
      
]
;   ]
   
end
   
   
be onePoint
      
ifelse != 99999999
      
[   if ok 
         
[
            
PenDown
            
setPC HSB kz 1 0.7
         
]
         
setXYZ kx kz y
         
PenUp
         
ok=true
      
][   ok=false
         
PenUp
      
]
   
end
   
ignore [
   
for [-360 360 360]
   
[   for [-1 1 2]
      
[   for [kz -100 100 dkz_s]
         
[   for [kx -100 100 dkx_s]
            
[   oneQuad
            
]            
         
]
      
]
   
]
]   
   
be oneQuad Ef
      
y1=E_solved_to_ky Ef kx kz d f
      
ifelse (abs dkx2) > 0.001 
      
and2 (abs (y1-oy)/dkx2) > dkx2min_s
      
and2 (abs y1-oy) < 100
      
[   dkx2=dkx2min_s*dkx2/(abs y1-oy)
      
][   dkx2=dkx_s
      
]
      
y2=E_solved_to_ky Ef kx+dkx2 kz d f
      
y3=E_solved_to_ky Ef kx+dkx2 kz+dkz_s d f
      
y4=E_solved_to_ky Ef kx kz+dkz_s d f
      
PenUp
      
setPC HSBA kz 1 0.7 1

ignore [
      
setXYZ kx kz y1
      
PolyStart
      
PenDown
      
setXYZ kx+dkx2 kz y2
      
setXYZ kx+dkx2 kz+dkz_s y3
      
setXYZ kx kz+dkz_s y4
      
setXYZ kx kz y1
      
PenUp
      
PolyEnd
]
;ignore[
      
if   y1 == 99999999 
      
and2 y2 == 99999999
      
and2 y3 == 99999999
      
and2 y4 == 99999999 [stop]
      
rim y1 == 99999999 
        
or2 y2 == 99999999
        
or2 y3 == 99999999
        
or2 y4 == 99999999
      
ifelse (mod round abs oy/100 2) == 0 
      
[   yh=d/a
      
][   yh=d/a+f*100
      
]
      
if y1 == 99999999 [y1=yh]
      
if y2 == 99999999 [y2=yh]
      
if y3 == 99999999 [y3=yh]
      
if y4 == 99999999 [y4=yh]
      
      
PolyStart
      
PenDown
      
setXYZ kx kz y1
      
setXYZ kx+dkx2 kz y2
      
setXYZ kx+dkx2 kz+dkz_s2 y3
      
setXYZ kx kz y1
      
PenUp
      
PolyEnd
      
      
PolyStart
      
PenDown
      
setXYZ kx kz y1
      
setXYZ kx+dkx2 kz+dkz_s2 y3
      
setXYZ kx kz+dkz_s2 y4
      
setXYZ kx kz y1
      
PenUp
      
PolyEnd
;]
      
oy=y1
   
end

;   setEye {0 -500 0}{0 0 0}{0 0 1}
;   redraw
   
drawEs
;   drawit

   
be drawEs
      
video=false
;      video=true
      
if video
      
[   (VideoStart "fermisurfaceE2.divx 25)
      
]
      
for [Ef -4 4 1/50]
      
[   clearScreen
         
draw Ef
         
updateGraph
         
if video [VideoFrame]
      
]
      
if video [VideoEnd]
   
end
end

to drawit
   
video=false
;   video=true
   
if video
   
[   (VideoStart "fermisurface.divx 25)
   
]
   
rotatescene3 video
   
if video [VideoEnd]
end

be rotatescene3 video [dphi 0.3][phi 0][dtheta 1][500]
   
local [eye r dr ddphi theta center upvector]
   
singleshot=Name? "framenr
   
if singleshot [phi=framenr*10]
   
eye=array 3
   
light=array 3
   
r=500
   
l=;200
   
dr=1.1
   
p={0 0 0}
   
dx={10 0 0}
   
theta=30
   
center={0 0 0}
   
upvector={0 1 0}
   
upv=upvector
   
ddphi=dphi/3
   
pr [leftright changes rotation speedup down set pitchESC exits]
   
dispatchMessages
   
forever
   
[   eye.1=r*(cos theta)*sin phi
      
eye.2=rsin theta
      
eye.3=r*(cos theta)*cos phi
      
upv.1=(cos 90-theta)*sin phi+180
      
upv.2sin 90-theta
      
upv.3=(cos 90-theta)*cos phi+180
      
setEye eye+center-p center-p upv
      
light.1=l*(cos theta)*sin phi
      
light.2=lsin theta
      
light.3=l*(cos theta)*cos phi
      
setLightPos light
      
redraw
      
if video [VideoFrame]
      
if singleshot [break]
      
phi=phi+dphi
      
if key? 
      
[   local [ch]
         
ch=readChar
         
ifelse ch==char 255
         
[   ch=readCharExt
            
ifElse (BitAnd MouseButtons 16)==16
            
[   dx=cross eye-center upvector
               
dx=(dx/Norm dx)*10
               
dy=cross dx eye-center
               
dy=(dy/Norm dy)*10
               
dz=eye-center
               
dz=(dz/Norm dz)*10
               
if ch==wxk_right [p=p+dx]
               
if ch==wxk_left  [p=p-dx]
               
if ch==wxk_up    [p=p+dy]
               
if ch==wxk_down  [p=p-dy]
               
if ch==wxk_prior [p=p+dz]
               
if ch==wxk_next  [p=p-dz]
            
][
               
if ch==wxk_right [dphi=dphi+ddphi]
               
if ch==wxk_left  [dphi=dphi-ddphi]
               
if ch==wxk_up    [theta=theta+dtheta]
               
if ch==wxk_down  [theta=theta-dtheta]
               
if ch==wxk_prior [r=r/dr]
               
if ch==wxk_next  [r=r*dr]
            
]
         
][
            
if ch==char 27 [stop]
            
if ch=="+ [r=r-dr]
            
if ch=="- [r=r+dr]
         
]
      
]
   
]
end