|
|
プログラムの後半部分です
---------------------
! ---------- 3D graphics
EXTERNAL SUB setRotateXYParameters(angleX,angleY,xCenter,yCenter,zCenter)
LET cosAx = COS(angleX)
LET sinAx = SIN(angleX)
LET cosAy = COS(angleY)
LET sinAy = SIN(angleY)
LET cx0 = xCenter
LET cy0 = yCenter
LET cz0 = zCenter
END SUB
EXTERNAL SUB plotLines3D(x,y,z,xShift,yShift) !shift*xRotateAx*yRotateAy*(shift^-1)
LET x1 = cosAy*(x-cx0)+sinAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
LET y1 = cosAx*(y-cy0)-sinAx*(z-cz0)
LET z1 =-sinAy*(x-cx0)+cosAy*(sinAx*(y-cy0)+cosAx*(z-cz0))
PLOT LINES: x1+cx0+xShift, y1+cy0+yShift; !z=z1+cz0
END SUB
! ---------- drawState
EXTERNAL SUB drawState(drawMode)
DECLARE EXTERNAL SUB setRotateXYParameters,drawState3D,draw3DPsix6,drawPsix6,drawRho
DECLARE EXTERNAL FUNCTION totalEnergy
SET DRAW MODE HIDDEN
CALL setRotateXYParameters(-PI/6,-PI/12,NNx/2,NNy/2,0)
CLEAR
SET TEXT HEIGHT 8
SET TEXT COLOR 1 ! black
IF drawMode=0 THEN
CALL drawRho(4,1.0,100,150) !(sc,zMag,xShift,yShift)
PLOT TEXT, AT 100, 150 :"electron density rho(x,y) and Veff(x,y)"
ELSEIF drawMode=1 THEN
CALL draw3DPsix6
ELSEIF drawMode=2 THEN
CALL drawPsix6
END IF
!--- caption
PLOT TEXT, AT 50,115 ,USING "iteration count =###### error =##.##########":iterCount,iterationError
PLOT TEXT, AT 50,100 ,USING "total energy =###.########## N of electron =##":totalEnergy,numberOfElectron
PLOT TEXT, AT 50, 85 ,USING "E0 =###.########## Occ =#.#####":sdEnergy(0),occupation(0)
PLOT TEXT, AT 50, 70 ,USING "E1 =###.########## Occ =#.#####":sdEnergy(1),occupation(1)
PLOT TEXT, AT 50, 55 ,USING "E2 =###.########## Occ =#.#####":sdEnergy(2),occupation(2)
PLOT TEXT, AT 50, 40 ,USING "E3 =###.########## Occ =#.#####":sdEnergy(3),occupation(3)
PLOT TEXT, AT 50, 25 ,USING "E4 =###.########## Occ =#.#####":sdEnergy(4),occupation(4)
PLOT TEXT, AT 50, 10 :"RS-DFT - Local Density Approximation 2D"
SET TEXT HEIGHT 10
PLOT TEXT, AT 50,470 :"[esc] exit [D] chage dispMode"
PLOT TEXT, AT 50,455 :"[1] ... [6] number of electron"
PLOT TEXT, AT 50,440 :"[7] hermonics k*x^2 [8] quantum well"
SET DRAW MODE EXPLICIT
END SUB
EXTERNAL SUB draw3DPsix6
DECLARE EXTERNAL SUB drawState3D
CALL drawState3D(0,2,0.5, 30,120) !(ist,sc,zMag,xShift,yShift)
CALL drawState3D(1,2,0.5,180,120)
CALL drawState3D(2,2,0.5,330,120)
CALL drawState3D(3,2,0.5, 30,270)
CALL drawState3D(4,2,0.5,180,270)
CALL drawState3D(5,2,0.5,330,270)
SET TEXT HEIGHT 10
PLOT TEXT, AT 30,140:"|0>"
PLOT TEXT, AT 180,140:"|1>"
PLOT TEXT, AT 330,140:"|2>"
PLOT TEXT, AT 30,290:"|3>"
PLOT TEXT, AT 180,290:"|4>"
PLOT TEXT, AT 330,290:"|5>"
END SUB
EXTERNAL SUB drawRho(sc,zMag,xShift,yShift)
DECLARE EXTERNAL SUB plotLines3D
FOR j=0 TO NNy-1 STEP 1
FOR i=0 TO NNx-1
LET psi = rho(i,j)*2000
IF psi>1 THEN
SET LINE COLOR 13
ELSE
SET LINE COLOR 3 ! potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),xShift,yShift) !(x,y,z)
NEXT i
PLOT LINES
NEXT j
FOR i=0 TO NNx-1 STEP 1
FOR j=0 TO NNy-1
LET psi = rho(i,j)*2000
IF psi>1 THEN
SET LINE COLOR 13
ELSE
SET LINE COLOR 3 ! potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),xShift,yShift) !(x,y,z)
NEXT j
PLOT LINES
NEXT i
END SUB
EXTERNAL SUB drawState3D(ist,sc,zMag,xShift,yShift)
DECLARE EXTERNAL SUB plotLines3D
FOR j=0 TO NNy-1 STEP 1
FOR i=0 TO NNx-1
LET psi = sdState(ist,i,j)*200
IF psi>1 THEN
SET LINE COLOR 4 ! red
ELSEIF psi<-1 THEN
SET LINE COLOR 2 ! blue
ELSE
SET LINE COLOR 3 ! potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),xShift,yShift) !(x,y,z)
NEXT i
PLOT LINES
NEXT j
FOR i=0 TO NNx-1 STEP 1
FOR j=0 TO NNy-1
LET psi = sdState(ist,i,j)*200
IF psi>1 THEN
SET LINE COLOR 4 ! red
ELSEIF psi<-1 THEN
SET LINE COLOR 2 ! blue
ELSE
SET LINE COLOR 3 ! potential:green
END IF
CALL plotLines3D(i*sc,j*sc,zMag*(psi + vv(i,j)),xShift,yShift) !(x,y,z)
NEXT j
PLOT LINES
NEXT i
END SUB
EXTERNAL SUB drawPsix6
DECLARE EXTERNAL SUB drawPsi
CALL drawPsi(0, 30,135) !(ist,xPos,yPos)
CALL drawPsi(1,180,135)
CALL drawPsi(2,330,135)
CALL drawPsi(3, 30,285)
CALL drawPsi(4,180,285)
CALL drawPsi(5,330,285)
SET TEXT HEIGHT 10
PLOT TEXT, AT 30,140:"|0>"
PLOT TEXT, AT 180,140:"|1>"
PLOT TEXT, AT 330,140:"|2>"
PLOT TEXT, AT 30,290:"|3>"
PLOT TEXT, AT 180,290:"|4>"
PLOT TEXT, AT 330,290:"|5>"
END SUB
EXTERNAL SUB drawPsi(ist,xPos,yPos) !drawMode=0
MAT psicol = vcol !--- psicol(,) <- Vcol(,) setted SUB initDraw
FOR i=0 TO NNx-1 !--- set psicol(,) for MAT PLOT CELLS
FOR j=0 TO NNy-1
LET psi = sdState(ist,i,j)
IF abs(psi)>0.002 THEN
IF psi>=0 THEN
LET col = psi*10
IF col>1 THEN LET col = 1
IF col<0 THEN LET col = 0
LET psicol(i,j) = 40+INT(col*50)
ELSE
LET col = -psi*10
IF col>1 THEN LET col = 1
IF col<0 THEN LET col = 0
LET psicol(i,j) = 100+INT(col*50)
END IF
END IF
NEXT j
NEXT i
MAT PLOT CELLS,IN xPos,yPos; xPos+128,yPos+128 :psicol
END SUB
END MODULE
|
|