楼主
请教一个程序，做毕业论文，各位高手不吝指教
PROGRAM ENGBND
REAL * 8H(33,33),A(33),E(33),EV(33,33),W(33)
REAL K,EL(21,3),KX(21),MSTAR(21,3)
DIMENSION WF(21,4),XX(21)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAK(20),NVAL,NPTS
CHARACTER SELTYN,EE,EH,ELEC(21,3)
LOGICAL LOG,LOG1,LOG2,LOG3
DATA EE,EH/'E','H'/
PI=3.1415926
NK=21
NX=21
OPEN(6,FILE='',STATUS='NEW')
OPEN(11,FILE='BND',STATUS='NEW')
OPEN(12,FILE='EMS',STATUS='NEW')
OPEN(13,FILE='PDF',STATUS='NEW')
WRITE(*,10)
WRITE(6,10)
10 FORMAT(' RUN THE PROGRAM TO CALCULATE',
/'FIRST THREE ENERGY LEVELS OF AN'/'ELE',
/CTRON FOR A GIVEN ERIODIC POTENTIAL.'/)
20 WRITE(*,30)
30 FORMAT('WHICH TYPE OF POTENTIALS DO YOU',
/'SELECT?'/' 1.RECTANGULAR;2.SAWTOOTH;3.'
/.'COSINE;4.HARMONIC;5.INTERPOLATED.'/2X,
/'INPUT THE NUMBER OF POTENTIALS(1-5):'\)
READ(*,40)IPTYPE
40 FORMAT(12)
GO TO(50,70,180,200,250)IPTYPE
50 WRITE(*,80)
WRITE(6,80)
60 FORMAT(/'THE RECTANGULAL POTENTIAL IS USED.')
GO TO 90
70 WRITE(*.80)
WRITE(6,80)
80 FORMAT(/’ THE SAWTOOTH POTENTIAL IS USED.’)
90 WRITE(*,100)
100 RMAT(‘ INPUT TH HEIGHT OF POTENTIAL:’
/ ,4X,’V0=’,\)
READ(*,*) V0
WRITE(6,110) V0
110 FORMAT(‘ THE HEIGHT OF POTENTIAL:’,4X,
/ ‘V0=’,F5.1)
IF(IPTYPE.NE.1) GO TO 140
WRITE(*,120)
120 ORMAT(‘ INPUT THE WIDTH OF RECTANGULAR:’
/ ,4X,’b=’,\)
READ(*,*) WIDTH
WTITE(6,130) WIDTH
130 FORMAT(‘THE WIDTH OF RECTANGULAR:’,
/ 4X,’b=’,F5.1)
140 WRITE(*,150)
150 FORMAT(‘ INPUT THE PERIOD OF POTENTIAL:’
/ ,5X,’a=’,\)
READ(*,*) PERIOD
WRITE(6,160) PERIOD
160 FORMAT(‘ THE PERIOD OF POTENTIAL:’,5X,
/ ‘a=’,F5.1)
IF(IPTYPE.NE.1.OR.PERIOD.GT.WIDTH) GOTO 380
WRITE(*,170)
WRITE(6,170)
170 FORMAT(‘PERIOD MUST BE GREATER THAN’
/ ,’WIDTH OF RECTANGLE!’)
GO TO 50
180 WRITE(*,190)
WRITE(6,190)
190 FORMAT(/’THE COSINE POTENTIAL IS USED,’,
/’IT PROPORTINAL TO 1-COS(2*PI*X/PERIOD).’)
GO TO 220
200 WRITE(*,210)
WRITE(*,210)
210 FORMAT(‘ INPUT CONSTANT OF’,
/’IS USED IT IS PROPORTIONALITY TO X**2.’)
220 WRITE(*,230)
230 FORMAT(‘ INPUT CONSTANT OF’,
/ ‘PROPORTIONALITY:V0=’/)
READ(*,*) V0
WRITE(6,240) V0
240 FORMAT(‘ VONSTANT O PROPORTIONALITY’
/ ,’:V0=’,F5.1)
WRITE(*,150)
READ(*,*) PERIOD
WRITE(6,160) PERIOD
I(IPTYPE.EQ.4) V0=V0*4.0/(PERIOD*PERIOD)
GO TO 380
250 WRITE(*,260)
WRITE(6,260)
260 FORMAT(/’ THE INTERPOLATED POTENTIAL’
/,’IS USED,WHICH IS GIVEN BY LINEAR’/
/’INTERPOLATION BETWEEN POINTS SPECIFIED’
/.’OVER HALF A PERIOD.’/’ THE PERIOD’
/.’IS TWICE THE FINAL X VALUE.’)
270 WRITE(*,280)
280 ORMAT(/’ INPUT THE NUMBER OF POINTS.’
/ ,11X,’NVAL=’./)
READ(*,40) NVAL
WRITE(6,290) NVAL
290 FORMAT(‘ THE NUMBER OF POINTS,’11X,
/ ‘NVAL=’,12)
IF(NVAL.GT.2.OR.NVAL.LT.20) GO TO 310
WRITE(*,300)
WRITE(6,300)
300 FORMAT(‘ THE NUMBER OF POINTS MUST BE’
/ ,’BETWEEN 2AND 201’)
GO TO 270
XVAL(1)=0.0
310 WRITE(*,350) 1,0.0,1
READ(*,*) FVAL(1)
DO 330 I=2,NVAL
WTITE(*,340) I,I
READ(*,) XVAL(I)
IF(XVAL(I).GT.XVAL(I-1)) GO TO 320
WRITE(*,370)
WRITE(6,370)
GO TO 310
320 WRITE(*,350) I,XVAL(I),I
READ(*,*) FVAL(I)
WRITE(6,360) I,XVAL(I),I,FVAL(I)
330 CONTINUE
340 FORMAT(‘ INPUT THE’,I2,’-TH VALUE’
/ ,’OF X,’,10X,’X(‘,I2,’)=’,\)
350 FORMAT(‘ INPUT THE POTENTIAL AT X(‘
/ ,I2,’)=’,F6.3,’,V(‘,I2,’)=’,\)
360 FORMAT(‘ THE POTENTIAL AT X(‘,I2,
/ ‘)=’,F6.3.’,V(‘,I2,’)=’,F6.3)
370 FORMAT(‘ EACH X VALUE MUST BE GREATER THAN LAST ONE!START AGAIN.’)
PERIOD=2.0*XVAL(NVAL)
WRITE(6,160) PERIOD
NPTS=100.0*PERIOD
380 CONTINUE
CALL FOURIR(33,DC,A)
WRITE(*,410)
WTITE(6,410)
410 FORMAT(/’THE FIRST THREE ENERGY LEVELS’,’FOR A SET OF WAVE VECTORS:’)
WRITE(11,420)NK
420 FORMAT(3X,I5)
WRITE(*,430)
WRITE(6,430)
430 FORMAT(/7X,’K’,9X,’E1’,8X,’E2’,8X,’E3’)
DO 550 L=1,NK
K=0.05*PI*FLOAT(L-1)/PERIOD
LOG=.FALSE.
IF(L.LE.10)NAR=10
IF(L.GT.10)NAR=11
440 CONTINUE
CALL ENERGY(K,E,EV,NAR,DC,A,H,W)
IF(LOG) GOTO 450
E1=E(1)
E2=E(2)
E3=E(3)
NAR=NAR+2
LOG=.TRUE.
GOTO 440
450 CONTINUE
F1=E(1)
F2=E(2)
F3=E(3)
LOG=.FALSE.
DE=F3-F1
IF(DE.EQ.0.0) DE=F3
NAR=NAR+8
LOG1=.TRUE.
LOG2=.TRUE.
LOG3=.TRUE.
ERES=ABS((F1-E1)/DE)
IF(ERES.GT.3.0E-3) LOG1=.FALSE.
ERES=ABS((F2-E2)/DE)
IF(ERES.GT.3.0E-3) LOG2=.FALSE.
ERES=ABS((F3-E3)/DE)
IF(ERES.GT.3.0E-3) LOG3=.FALSE.
IF(LOG1.AND.LOG2.AND.LOG3) GOTO 530
IF(NAR.LE.31)GOTO 440
IF(.NOT.LOG1.AND.LOG2.AND.LOG3) WRITE(*,460)
WRITE(6,800)
800 FORMAT(6X,’NOTE: E---ELECTRON,’,’H---HOLE.’)
WRITE(12,420)NK
DO 820 L=1,NK
DO 810 IN=1,3
810 IF(MSTAR(L,IN).GT.50.0) MSTAR(L,IN)=50.0
WRITE(12,560) KX(L),(MSTAR(L,IN),IN=1,3)
820 CONTINUE
CLOSE(12)
830 WRITE(*,840)
840 FORMAT(/’ DO YOU WANT TO CALCULATE THE’
/ ,’PROBABILITY-DENSITY FUNTIONS’/’FOR’
/ ,’THE FIRST THREE ENERGY STATESLY/N)’\)
READ(*,580) SELTYN
IF(SELTYN.NE.’Y’.AND.SELTYN.NE.’y’) GOTO 980
910 WRITE(*,920)
920 FORMAT(‘ INPUT THE VALUE OF K AS AREAL’
/ ,’NUMBER,K+’,\)
READ(*,*) K
WRITE(*,930) K
WRITE(6,930) K
930 FORMAT(/’ THE PROBABILITY-DENSITY FUNTIONS’
/ ,’OF THE FIRST THREE ENERGY STATES:’,/,
/ ‘(K=’,F6.4,’)’)
CALL ENERGY(K,E,EV,15,DC,A,H,W)
DO 950 LX=1,NX
XX(LX)=0.05*FLOAT(LX-11)*PERIOD
X=XX(LX)
WF(LX,4)=V(X)
DO 940 IN=1,3
WF(LX,IN)=PROB(IN,EV,15,X)
940 CONTINUE
950 CONTINUE
WRITE(13,420) NX
WRITE(*,960)
WRITE(6,960)
960 FORMAT(/7X,’X’,9X,’PD1’,7X,’PD2’,7X,’PD3’,7X,’V’)
DO 970 LX=1,NX
WRITE(*,560) XX(LX),(WF(LX,IN),IN=1,4)
WTITE(6,560) XX(LX),( WF(LX,IN),IN=1,4)
WRITE(13,560) XX(LX),(WF(LX,IN),IN=1,4)
970 CONTINUE
CLOSE(13)
980 STOP
END
SUBROUTINE FOURIR(NFOUR,DC,A)
REAL *8 A(NFOUR)
REAL K
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(200,XVAL(20),NVAL,NPTS
GO TO(10,30,50,70,90) IPTYPE
10 DC=WIDTH*V0/PERIOD
DO 20 M=1,NFOUR
FM=FLOAT(M)
A(M)=2*COS(FM*PI)*V0*SIN(FM*PI*WIDTH/PERIOD)/(M*PI)
20 CONTINUE
RETURN
30 DC=V0/2.0
DO 40 M=1,NFOUR
FM=FLOAT(M)
A(M)=2.0*V0*(COS(FM*PI)-1.0)/(FM*PI)**2
40 CONTINUE
RETURN
50 DC=V0
A(1)=-V0
DO 60 M=2,NFOUR
A(M)=0.0
60 CONTINUE
RETURN
70 DC=V0*PERIOD**2/12.0
DO 80 M=1,NFOUR
FM=FLOAT(M)
A(M)=V0*PERIOD**2*COS(FM*PI)/(FM*PI)**2
80 CONTINUE
RETURN
90 K=2.0*PI/PERIOD
FD=FLOAT(NPTS)
STEP=PERIOS/FD
DC=V(0.0)
100 CONTINUE
DO 120 I=1,NPTS-1
X=FLOAT(I)*STEP
VX=V(X)
DC=DC+VX
DO 110 M=1,NFOUR
FM=FLOAT(M)
A(M)=A(M)+COS(FM*K*X)*VX
110 CONTINUE
120 CONTINUE
DC=DC/FD
DO 130 M=1,NFOUR
A(M)=A(M)*2.0/FD
130 CONTINUE
RETURN
END
SUBROUTINE ENERGY(K,E,EV,N,DC,A,H,W)
REAL *8 E(N),H(N,N),A(N),EV(N,N),W(N)
REAL K
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAL(20),NVAL,NPTS
N1=N/2
IF(2*N1.EQ.N) GOTO 20
N2=(N+1)/2
DO 10 I=1,N
H(I,I)=(K-FLOAT(N2-I)*2.0*PI/PERIOD)**2+DC
10 CONTINUE
GOTO 40
20 DO 30 I=1,N
H(I,I)=( K-FLOAT(N1-I+1)*2.0*PI/PERIOD)**2+DC
30 CONTINUE
40 DO 50 I=1,N+1
DO 50 J=I+1,N
JMI=J-I
H(I,J)=A(JMI)/2.0
H(J,I)=H(I,J)
50 CONTINUE
CALL EIGRAR(H,E,EV,W,N)
RETURN
END
FUNCTION V(X)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAL(20),NVAL,NPTS
IF(X.LT.0.0) X=-X
P2=PERIOD/2.0
R=AMOD(X,PERIOD)
IF(R.GT.P2) R=PERIOD-R
GO TO (10,20,30,40,50) IPTYPE
10 B=(PERIOD-WIDTH)/2.0
IF(R.LT.B) V=0.0
IF(R.EQ.B) V=V0/2.0
IF(R.GT.B) V=V0
RETURN
20 V=R*V0/P2
RETURN
30 Q=2.0*PI/PERIOD
V=V0*(1.0-COS(Q*R))
RETURN
40 V=V0*R**2
RETURN
50 NV1=NVAL-1
DO 60 I=1,NV1
IF(R.GE.XVAL(I).AND.R.LT.XVAL(I+1)) J=1
60 CONTINUE
GRAD=(FVAL(J+1)-FVAL(J))/(XVAL(J+1)-XVAL(J))
V=VAL(J)+GRAD*(R-XVAL(J))
RETURN
END
FUNCTION PROB(L,EV,N,X)
REAL *8 EV(N,N)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,VAL(20),XVAL(20),NVAL,NPTS
S=0.0
DO 20 I=1,N-1
DO 10 J=I+1,N
F=2.0*PI*FLOAT(J-I)/PERIOD
S=S+EV(I,L)*EV(J,L)*CLOS(F*X)
10 CONTINUE
20 CONTINUE
PROB=1.0+2.0*S
RETURN
END
SUBROUTINE EIGRAR(H,EA,EV,W,N)
REAL *8 H(N,N),EA(N),W(N)
CALL MATEIG(H,EV,N)
DO 10 I=1,N
EA(I)=H(I,I)
10 CONTINUE
DO 40 JJ=1,N-1
DO 30 KK=JJ+1,N
IF((EA(JJ)-EA(KK)).LE.0.0) GOTO 30
W(JJ)=EA(JJ)
EA(JJ)=EA(KK)
DO 20 LN=1,N
W(LN)=EV(LN,JJ)
EV(LN,JJ)=EV(LN,KK)
EV(LN,KK)=W(LN)
20 CONTINUE
30 CONTINUE
40 CONTINUE
RETURN
END
SUBROUTINE MATEIG(H,EV,N)
REAL *8 H(N,N),EV(N,N),ANORM,ETC,CONSTF,B,C
REAL *8 UA,VA,WA,SNN,SNN2,CSN,CSN2,APP,AQQ
EPS=0.1D-5
CONSTF=FLOAT(N)
IN=0
DO 30 J=1,N
DO 20 J=1,N
I((I-J).EQ.0) GOTO 10
EV(I,J)=0.0
GOTO 20
10 EV(I,J)=1.0
20 CONTINUE
30 CONTINUE
ANORM=0.0
DO 50 J=1,N
DO 40 I=1,N
IF((J-I).EQ.0) GOTO 40
ANORM=ANORM+H(I,J)*H(I,J)
40 CONTINUE
50 CONTINUE
ANORM=DSQRT(ANORM)
FNORM=ANORM*EPS
60 ANORM=ANORM/CONSTF
70 DO 110 IQ=2,N
DO 110 IP=1,IQ-1
IF((DABS(H(IP,IQ))-ANORM).LE.00) GO TO 110
IN=1
UA=-H(IP,IQ)
VA=(H(IP,IP)-H(IQ,IQ))/2.0
WA=UA/(DSQRT(UA*UA+VA*VA))
IF(VA.LE.0.0) WA=-WA
SNN=WA/(DSQRT(2.0*(1.0+DSQRT(1.0-WA*WA))))
SNN2=SNN*SNN
CSN=DSQRT(1.0-SNN2)
DO 90 I=1,N
IF ((I-IP).EQ.0) GO TO 80
IF ((I-IQ).EQ.0) GO TO 80
B=H(I,IP)*CSN-H(I,IQ)*SNN
H(I,IQ)=H(I,IP)*SNN+H(I,IQ)*CSN
H(I,IP)=B
80 CONTINUE
C=EV(I,IP)*CSN-EV(I,IQ)*SNN
EV(I,IQ)=EV(I,IQ)*CSN+EV(I,IP)*SNN
EV(I,IP)=C
90 CONTINUE
CSN2=CSN*CSN
STC=SNN*CSN
APP=H(IP,IP)*CSN2+H(IQ,IQ)*SNN2-2.0*H(IP,IQ)*STC
AQQ=H(IP,IP)*SNN2+H(IQ,IQ)8CSN2+2.0*H(IP,IQ)*STC
H(IP,IQ)=(H(IP,IP)-H(IQ,IQ))*STC+H(IP,IQ)*(CSN2-SNN2)
H(IQ,IP)=H(IP,IQ)
H(IP,IP)=APP
H(IQ,IQ)=AQQ
DO 100 I=1,N
H(IP,I)=H(I,IP)
100 H(IQ,I)=H(I,IQ)
110 CONTINUE
IF((IN-1).NE.0) GO TO 120
IN=0
GO TO 70
120 IF((ANORM-FNORM).GT.0.0) GO TO 60
RETURN
END
PROGRAM ENGBND
REAL * 8H(33,33),A(33),E(33),EV(33,33),W(33)
REAL K,EL(21,3),KX(21),MSTAR(21,3)
DIMENSION WF(21,4),XX(21)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAK(20),NVAL,NPTS
CHARACTER SELTYN,EE,EH,ELEC(21,3)
LOGICAL LOG,LOG1,LOG2,LOG3
DATA EE,EH/'E','H'/
PI=3.1415926
NK=21
NX=21
OPEN(6,FILE='',STATUS='NEW')
OPEN(11,FILE='BND',STATUS='NEW')
OPEN(12,FILE='EMS',STATUS='NEW')
OPEN(13,FILE='PDF',STATUS='NEW')
WRITE(*,10)
WRITE(6,10)
10 FORMAT(' RUN THE PROGRAM TO CALCULATE',
/'FIRST THREE ENERGY LEVELS OF AN'/'ELE',
/CTRON FOR A GIVEN ERIODIC POTENTIAL.'/)
20 WRITE(*,30)
30 FORMAT('WHICH TYPE OF POTENTIALS DO YOU',
/'SELECT?'/' 1.RECTANGULAR;2.SAWTOOTH;3.'
/.'COSINE;4.HARMONIC;5.INTERPOLATED.'/2X,
/'INPUT THE NUMBER OF POTENTIALS(1-5):'\)
READ(*,40)IPTYPE
40 FORMAT(12)
GO TO(50,70,180,200,250)IPTYPE
50 WRITE(*,80)
WRITE(6,80)
60 FORMAT(/'THE RECTANGULAL POTENTIAL IS USED.')
GO TO 90
70 WRITE(*.80)
WRITE(6,80)
80 FORMAT(/’ THE SAWTOOTH POTENTIAL IS USED.’)
90 WRITE(*,100)
100 RMAT(‘ INPUT TH HEIGHT OF POTENTIAL:’
/ ,4X,’V0=’,\)
READ(*,*) V0
WRITE(6,110) V0
110 FORMAT(‘ THE HEIGHT OF POTENTIAL:’,4X,
/ ‘V0=’,F5.1)
IF(IPTYPE.NE.1) GO TO 140
WRITE(*,120)
120 ORMAT(‘ INPUT THE WIDTH OF RECTANGULAR:’
/ ,4X,’b=’,\)
READ(*,*) WIDTH
WTITE(6,130) WIDTH
130 FORMAT(‘THE WIDTH OF RECTANGULAR:’,
/ 4X,’b=’,F5.1)
140 WRITE(*,150)
150 FORMAT(‘ INPUT THE PERIOD OF POTENTIAL:’
/ ,5X,’a=’,\)
READ(*,*) PERIOD
WRITE(6,160) PERIOD
160 FORMAT(‘ THE PERIOD OF POTENTIAL:’,5X,
/ ‘a=’,F5.1)
IF(IPTYPE.NE.1.OR.PERIOD.GT.WIDTH) GOTO 380
WRITE(*,170)
WRITE(6,170)
170 FORMAT(‘PERIOD MUST BE GREATER THAN’
/ ,’WIDTH OF RECTANGLE!’)
GO TO 50
180 WRITE(*,190)
WRITE(6,190)
190 FORMAT(/’THE COSINE POTENTIAL IS USED,’,
/’IT PROPORTINAL TO 1-COS(2*PI*X/PERIOD).’)
GO TO 220
200 WRITE(*,210)
WRITE(*,210)
210 FORMAT(‘ INPUT CONSTANT OF’,
/’IS USED IT IS PROPORTIONALITY TO X**2.’)
220 WRITE(*,230)
230 FORMAT(‘ INPUT CONSTANT OF’,
/ ‘PROPORTIONALITY:V0=’/)
READ(*,*) V0
WRITE(6,240) V0
240 FORMAT(‘ VONSTANT O PROPORTIONALITY’
/ ,’:V0=’,F5.1)
WRITE(*,150)
READ(*,*) PERIOD
WRITE(6,160) PERIOD
I(IPTYPE.EQ.4) V0=V0*4.0/(PERIOD*PERIOD)
GO TO 380
250 WRITE(*,260)
WRITE(6,260)
260 FORMAT(/’ THE INTERPOLATED POTENTIAL’
/,’IS USED,WHICH IS GIVEN BY LINEAR’/
/’INTERPOLATION BETWEEN POINTS SPECIFIED’
/.’OVER HALF A PERIOD.’/’ THE PERIOD’
/.’IS TWICE THE FINAL X VALUE.’)
270 WRITE(*,280)
280 ORMAT(/’ INPUT THE NUMBER OF POINTS.’
/ ,11X,’NVAL=’./)
READ(*,40) NVAL
WRITE(6,290) NVAL
290 FORMAT(‘ THE NUMBER OF POINTS,’11X,
/ ‘NVAL=’,12)
IF(NVAL.GT.2.OR.NVAL.LT.20) GO TO 310
WRITE(*,300)
WRITE(6,300)
300 FORMAT(‘ THE NUMBER OF POINTS MUST BE’
/ ,’BETWEEN 2AND 201’)
GO TO 270
XVAL(1)=0.0
310 WRITE(*,350) 1,0.0,1
READ(*,*) FVAL(1)
DO 330 I=2,NVAL
WTITE(*,340) I,I
READ(*,) XVAL(I)
IF(XVAL(I).GT.XVAL(I-1)) GO TO 320
WRITE(*,370)
WRITE(6,370)
GO TO 310
320 WRITE(*,350) I,XVAL(I),I
READ(*,*) FVAL(I)
WRITE(6,360) I,XVAL(I),I,FVAL(I)
330 CONTINUE
340 FORMAT(‘ INPUT THE’,I2,’-TH VALUE’
/ ,’OF X,’,10X,’X(‘,I2,’)=’,\)
350 FORMAT(‘ INPUT THE POTENTIAL AT X(‘
/ ,I2,’)=’,F6.3,’,V(‘,I2,’)=’,\)
360 FORMAT(‘ THE POTENTIAL AT X(‘,I2,
/ ‘)=’,F6.3.’,V(‘,I2,’)=’,F6.3)
370 FORMAT(‘ EACH X VALUE MUST BE GREATER THAN LAST ONE!START AGAIN.’)
PERIOD=2.0*XVAL(NVAL)
WRITE(6,160) PERIOD
NPTS=100.0*PERIOD
380 CONTINUE
CALL FOURIR(33,DC,A)
WRITE(*,410)
WTITE(6,410)
410 FORMAT(/’THE FIRST THREE ENERGY LEVELS’,’FOR A SET OF WAVE VECTORS:’)
WRITE(11,420)NK
420 FORMAT(3X,I5)
WRITE(*,430)
WRITE(6,430)
430 FORMAT(/7X,’K’,9X,’E1’,8X,’E2’,8X,’E3’)
DO 550 L=1,NK
K=0.05*PI*FLOAT(L-1)/PERIOD
LOG=.FALSE.
IF(L.LE.10)NAR=10
IF(L.GT.10)NAR=11
440 CONTINUE
CALL ENERGY(K,E,EV,NAR,DC,A,H,W)
IF(LOG) GOTO 450
E1=E(1)
E2=E(2)
E3=E(3)
NAR=NAR+2
LOG=.TRUE.
GOTO 440
450 CONTINUE
F1=E(1)
F2=E(2)
F3=E(3)
LOG=.FALSE.
DE=F3-F1
IF(DE.EQ.0.0) DE=F3
NAR=NAR+8
LOG1=.TRUE.
LOG2=.TRUE.
LOG3=.TRUE.
ERES=ABS((F1-E1)/DE)
IF(ERES.GT.3.0E-3) LOG1=.FALSE.
ERES=ABS((F2-E2)/DE)
IF(ERES.GT.3.0E-3) LOG2=.FALSE.
ERES=ABS((F3-E3)/DE)
IF(ERES.GT.3.0E-3) LOG3=.FALSE.
IF(LOG1.AND.LOG2.AND.LOG3) GOTO 530
IF(NAR.LE.31)GOTO 440
IF(.NOT.LOG1.AND.LOG2.AND.LOG3) WRITE(*,460)
WRITE(6,800)
800 FORMAT(6X,’NOTE: E---ELECTRON,’,’H---HOLE.’)
WRITE(12,420)NK
DO 820 L=1,NK
DO 810 IN=1,3
810 IF(MSTAR(L,IN).GT.50.0) MSTAR(L,IN)=50.0
WRITE(12,560) KX(L),(MSTAR(L,IN),IN=1,3)
820 CONTINUE
CLOSE(12)
830 WRITE(*,840)
840 FORMAT(/’ DO YOU WANT TO CALCULATE THE’
/ ,’PROBABILITY-DENSITY FUNTIONS’/’FOR’
/ ,’THE FIRST THREE ENERGY STATESLY/N)’\)
READ(*,580) SELTYN
IF(SELTYN.NE.’Y’.AND.SELTYN.NE.’y’) GOTO 980
910 WRITE(*,920)
920 FORMAT(‘ INPUT THE VALUE OF K AS AREAL’
/ ,’NUMBER,K+’,\)
READ(*,*) K
WRITE(*,930) K
WRITE(6,930) K
930 FORMAT(/’ THE PROBABILITY-DENSITY FUNTIONS’
/ ,’OF THE FIRST THREE ENERGY STATES:’,/,
/ ‘(K=’,F6.4,’)’)
CALL ENERGY(K,E,EV,15,DC,A,H,W)
DO 950 LX=1,NX
XX(LX)=0.05*FLOAT(LX-11)*PERIOD
X=XX(LX)
WF(LX,4)=V(X)
DO 940 IN=1,3
WF(LX,IN)=PROB(IN,EV,15,X)
940 CONTINUE
950 CONTINUE
WRITE(13,420) NX
WRITE(*,960)
WRITE(6,960)
960 FORMAT(/7X,’X’,9X,’PD1’,7X,’PD2’,7X,’PD3’,7X,’V’)
DO 970 LX=1,NX
WRITE(*,560) XX(LX),(WF(LX,IN),IN=1,4)
WTITE(6,560) XX(LX),( WF(LX,IN),IN=1,4)
WRITE(13,560) XX(LX),(WF(LX,IN),IN=1,4)
970 CONTINUE
CLOSE(13)
980 STOP
END
SUBROUTINE FOURIR(NFOUR,DC,A)
REAL *8 A(NFOUR)
REAL K
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(200,XVAL(20),NVAL,NPTS
GO TO(10,30,50,70,90) IPTYPE
10 DC=WIDTH*V0/PERIOD
DO 20 M=1,NFOUR
FM=FLOAT(M)
A(M)=2*COS(FM*PI)*V0*SIN(FM*PI*WIDTH/PERIOD)/(M*PI)
20 CONTINUE
RETURN
30 DC=V0/2.0
DO 40 M=1,NFOUR
FM=FLOAT(M)
A(M)=2.0*V0*(COS(FM*PI)-1.0)/(FM*PI)**2
40 CONTINUE
RETURN
50 DC=V0
A(1)=-V0
DO 60 M=2,NFOUR
A(M)=0.0
60 CONTINUE
RETURN
70 DC=V0*PERIOD**2/12.0
DO 80 M=1,NFOUR
FM=FLOAT(M)
A(M)=V0*PERIOD**2*COS(FM*PI)/(FM*PI)**2
80 CONTINUE
RETURN
90 K=2.0*PI/PERIOD
FD=FLOAT(NPTS)
STEP=PERIOS/FD
DC=V(0.0)
100 CONTINUE
DO 120 I=1,NPTS-1
X=FLOAT(I)*STEP
VX=V(X)
DC=DC+VX
DO 110 M=1,NFOUR
FM=FLOAT(M)
A(M)=A(M)+COS(FM*K*X)*VX
110 CONTINUE
120 CONTINUE
DC=DC/FD
DO 130 M=1,NFOUR
A(M)=A(M)*2.0/FD
130 CONTINUE
RETURN
END
SUBROUTINE ENERGY(K,E,EV,N,DC,A,H,W)
REAL *8 E(N),H(N,N),A(N),EV(N,N),W(N)
REAL K
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAL(20),NVAL,NPTS
N1=N/2
IF(2*N1.EQ.N) GOTO 20
N2=(N+1)/2
DO 10 I=1,N
H(I,I)=(K-FLOAT(N2-I)*2.0*PI/PERIOD)**2+DC
10 CONTINUE
GOTO 40
20 DO 30 I=1,N
H(I,I)=( K-FLOAT(N1-I+1)*2.0*PI/PERIOD)**2+DC
30 CONTINUE
40 DO 50 I=1,N+1
DO 50 J=I+1,N
JMI=J-I
H(I,J)=A(JMI)/2.0
H(J,I)=H(I,J)
50 CONTINUE
CALL EIGRAR(H,E,EV,W,N)
RETURN
END
FUNCTION V(X)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,FVAL(20),XVAL(20),NVAL,NPTS
IF(X.LT.0.0) X=-X
P2=PERIOD/2.0
R=AMOD(X,PERIOD)
IF(R.GT.P2) R=PERIOD-R
GO TO (10,20,30,40,50) IPTYPE
10 B=(PERIOD-WIDTH)/2.0
IF(R.LT.B) V=0.0
IF(R.EQ.B) V=V0/2.0
IF(R.GT.B) V=V0
RETURN
20 V=R*V0/P2
RETURN
30 Q=2.0*PI/PERIOD
V=V0*(1.0-COS(Q*R))
RETURN
40 V=V0*R**2
RETURN
50 NV1=NVAL-1
DO 60 I=1,NV1
IF(R.GE.XVAL(I).AND.R.LT.XVAL(I+1)) J=1
60 CONTINUE
GRAD=(FVAL(J+1)-FVAL(J))/(XVAL(J+1)-XVAL(J))
V=VAL(J)+GRAD*(R-XVAL(J))
RETURN
END
FUNCTION PROB(L,EV,N,X)
REAL *8 EV(N,N)
COMMON V0,WIDTH,PERIOD,PI,IPTYPE,VAL(20),XVAL(20),NVAL,NPTS
S=0.0
DO 20 I=1,N-1
DO 10 J=I+1,N
F=2.0*PI*FLOAT(J-I)/PERIOD
S=S+EV(I,L)*EV(J,L)*CLOS(F*X)
10 CONTINUE
20 CONTINUE
PROB=1.0+2.0*S
RETURN
END
SUBROUTINE EIGRAR(H,EA,EV,W,N)
REAL *8 H(N,N),EA(N),W(N)
CALL MATEIG(H,EV,N)
DO 10 I=1,N
EA(I)=H(I,I)
10 CONTINUE
DO 40 JJ=1,N-1
DO 30 KK=JJ+1,N
IF((EA(JJ)-EA(KK)).LE.0.0) GOTO 30
W(JJ)=EA(JJ)
EA(JJ)=EA(KK)
DO 20 LN=1,N
W(LN)=EV(LN,JJ)
EV(LN,JJ)=EV(LN,KK)
EV(LN,KK)=W(LN)
20 CONTINUE
30 CONTINUE
40 CONTINUE
RETURN
END
SUBROUTINE MATEIG(H,EV,N)
REAL *8 H(N,N),EV(N,N),ANORM,ETC,CONSTF,B,C
REAL *8 UA,VA,WA,SNN,SNN2,CSN,CSN2,APP,AQQ
EPS=0.1D-5
CONSTF=FLOAT(N)
IN=0
DO 30 J=1,N
DO 20 J=1,N
I((I-J).EQ.0) GOTO 10
EV(I,J)=0.0
GOTO 20
10 EV(I,J)=1.0
20 CONTINUE
30 CONTINUE
ANORM=0.0
DO 50 J=1,N
DO 40 I=1,N
IF((J-I).EQ.0) GOTO 40
ANORM=ANORM+H(I,J)*H(I,J)
40 CONTINUE
50 CONTINUE
ANORM=DSQRT(ANORM)
FNORM=ANORM*EPS
60 ANORM=ANORM/CONSTF
70 DO 110 IQ=2,N
DO 110 IP=1,IQ-1
IF((DABS(H(IP,IQ))-ANORM).LE.00) GO TO 110
IN=1
UA=-H(IP,IQ)
VA=(H(IP,IP)-H(IQ,IQ))/2.0
WA=UA/(DSQRT(UA*UA+VA*VA))
IF(VA.LE.0.0) WA=-WA
SNN=WA/(DSQRT(2.0*(1.0+DSQRT(1.0-WA*WA))))
SNN2=SNN*SNN
CSN=DSQRT(1.0-SNN2)
DO 90 I=1,N
IF ((I-IP).EQ.0) GO TO 80
IF ((I-IQ).EQ.0) GO TO 80
B=H(I,IP)*CSN-H(I,IQ)*SNN
H(I,IQ)=H(I,IP)*SNN+H(I,IQ)*CSN
H(I,IP)=B
80 CONTINUE
C=EV(I,IP)*CSN-EV(I,IQ)*SNN
EV(I,IQ)=EV(I,IQ)*CSN+EV(I,IP)*SNN
EV(I,IP)=C
90 CONTINUE
CSN2=CSN*CSN
STC=SNN*CSN
APP=H(IP,IP)*CSN2+H(IQ,IQ)*SNN2-2.0*H(IP,IQ)*STC
AQQ=H(IP,IP)*SNN2+H(IQ,IQ)8CSN2+2.0*H(IP,IQ)*STC
H(IP,IQ)=(H(IP,IP)-H(IQ,IQ))*STC+H(IP,IQ)*(CSN2-SNN2)
H(IQ,IP)=H(IP,IQ)
H(IP,IP)=APP
H(IQ,IQ)=AQQ
DO 100 I=1,N
H(IP,I)=H(I,IP)
100 H(IQ,I)=H(I,IQ)
110 CONTINUE
IF((IN-1).NE.0) GO TO 120
IN=0
GO TO 70
120 IF((ANORM-FNORM).GT.0.0) GO TO 60
RETURN
END
