DIMENSION X(1100),P(1100),H(1100),W(2200)
COMMON /COMAK/AK(0:1100)
DATA NW,PAI/2200,3.14159265/
OPEN(8,FILE='OUT.DAT',STATUS='UNKNOWN')
N=129
CALL SUBAK(N)
DX=3./(N-1)
DO I=1,N
X(I)=-1.5+(I-1)*DX
P(I)=0.0
IF(ABS(X(I)).LE.1.0)P(I)=SQRT(1.-X(I)*X(I))
ENDDO
K=3
CALL DISP(N,NW,K,DX,P,W)
WX=0.5*PAI*DX*ALOG(DX)
DO I=1,N
W(I)=WX
DO J=1,N
IJ=IABS(I-J)
W(I)=W(I)+AK(IJ)*P(J)*DX
ENDDO
ENDDO
DO 30 I=1,N
H(I)=1.24+0.5*X(I)*X(I)-W(I)/PAI
30 CONTINUE
DO I=1,N
WRITE(8,40)X(I),P(I),H(I)
ENDDO
40 FORMAT(1X,6(E12.6,1X))
STOP
END
SUBROUTINE DISP(N,NW,KMAX,DX,P1,W)
DIMENSION P1(N),W(NW),P(2200),AK1(0:50),AK2(0:50)
COMMON /COMAK/AK(0:1100)
DATA NMAX,KMIN/2200,1/
N2=N
M=3+2*ALOG(FLOAT(N))
K1=N+KMAX
DO 10 I=1,N
10 P(K1+I)=P1(I)
DO 20 KK=KMIN,KMAX-1
K=KMAX+KMIN-KK
N1=(N2+1)/2
CALL DOWNP(NMAX,N1,N2,K,P)
20 N2=N1
DX1=DX*2**(KMAX-KMIN)
CALL WI(NMAX,N1,KMIN,KMAX,DX,DX1,P,W)
DO 30 K=KMIN+1,KMAX
N2=2*N1-1
DX1=DX1/2.
CALL AKCO(M+5,KMAX,K,AK1)
CALL AKIN(M+6,AK1,AK2)
CALL WCOS(NMAX,N1,N2,K,W)
CALL CORR(NMAX,N2,K,M,1,DX1,P,W,AK1)
CALL WINT(NMAX,N2,K,W)
CALL CORR(NMAX,N2,K,M,2,DX1,P,W,AK2)
30 N1=N2
DO 40 I=1,N
40 W(I)=W(K1+I)
RETURN
END
SUBROUTINE DOWNP(NMAX,N1,N2,K,P)
DIMENSION P(NMAX)
K1=N1+K-1
K2=N2+K-1
DO 10 I=3,N1-2
I2=2*I+K2
10 P(K1+I)=(16.*P(I2)+9.*(P(I2-1)+P(I2+1))-(P(I2-3)+P(I2+3)))/32.
P(K1+2)=0.25*(P(K2+3)+P(K2+5))+0.5*P(K2+4)
P(K1+N1-1)=0.25*(P(K2+N2-2)+P(K2+N2))+0.5*P(K2+N2-1)
RETURN
END
SUBROUTINE WCOS(NMAX,N1,N2,K,W)
DIMENSION W(NMAX)
K1=N1+K-1
K2=N2+K
DO 10 I=1,N1
II=2*I-1
10 W(K2+II)=W(K1+I)
RETURN
END
SUBROUTINE WINT(NMAX,N,K,W)
DIMENSION W(NMAX)
K2=N+K
DO 10 I=4,N-3,2
II=K2+I
10 W(II)=(9.*(W(II-1)+W(II+1))-(W(II-3)+W(II+3)))/16.
I1=K2+2
I2=K2+N-1
W(I1)=0.5*(W(I1-1)+W(I1+1))
W(I2)=0.5*(W(I2-1)+W(I2+1))
RETURN
END
SUBROUTINE CORR(NMAX,N,K,M,I1,DX,P,W,AK)
DIMENSION P(NMAX),W(NMAX),AK(0:M)
K1=N+K
IF(I1.EQ.2)GOTO 20
DO 10 I=1,N,2
II=K1+I
J1=MAX0(1,I-M)
J2=MIN0(N,I+M)
DO 10 J=J1,J2
IJ=IABS(I-J)
10 W(II)=W(II)+AK(IJ)*DX*P(K1+J)
RETURN
20 DO 30 I=2,N,2
II=K1+I
J1=MAX0(1,I-M)
J2=MIN0(N,I+M)
DO 30 J=J1,J2
IJ=IABS(I-J)
30 W(II)=W(II)+AK(IJ)*DX*P(K1+J)
RETURN
END
SUBROUTINE WI(NMAX,N,KMIN,KMAX,DX,DX1,P,W)
DIMENSION P(NMAX),W(NMAX)
COMMON /COMAK/AK(0:1100)
K1=N+1
K=2**(KMAX-KMIN)
C=ALOG(DX)
DO 10 I=1,N
II=K1+I
W(II)=0.0
DO 10 J=1,N
IJ=K*IABS(I-J)
10 W(II)=W(II)+(AK(IJ)+C)*DX1*P(K1+J)
RETURN
END
SUBROUTINE AKCO(KA,KMAX,K,AK1)
DIMENSION AK1(0:KA)
COMMON /COMAK/AK(0:1100)
J=2**(KMAX-K)
DO 10 I=0,KA
II=J*I
10 AK1(I)=AK(II)
RETURN
END
SUBROUTINE AKIN(KA,AK1,AK2)
DIMENSION AK1(KA),AK2(KA)
DO 10 I=4,KA-3
10 AK2(I)=(9.*(AK1(I-1)+AK1(I+1))-(AK1(I-3)+AK1(I+3)))/16.
AK2(1)=(9.*AK1(2)-AK1(4))/8.
AK2(2)=(9.*(AK1(1)+AK1(3))-(AK1(3)+AK1(5)))/16.
AK2(3)=(9.*(AK1(2)+AK1(4))-(AK1(2)+AK1(6)))/16.
DO 20 I=1,KA
20 AK2(I)=AK1(I)-AK2(I)
DO 30 I=1,KA-1,2
I1=I+1
AK1(I)=0.0
30 AK1(I1)=AK2(I1)
RETURN
END
SUBROUTINE SUBAK(MM)
COMMON /COMAK/AK(0:1100)
DO 10 I=0,MM
10 AK(I)=(I+0.5)*(ALOG(ABS(I+0.5))-1.)-(I-0.5)*(ALOG(ABS(I-0.5))-1.)
RETURN
END |