Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- C Copy of INTEGRAL FORTRAN WRITTEN BY Jay Rasaiah
- C MODIFIED INPUT TO RUN ON IMAC by E.LOVEJOY June 18 2024, Rev 8.01
- C
- C Compiled with:
- C gfortran -g -fcheck=all -std=legacy -o integral integral.f
- C -fcheck=all enables array bound checks at run time
- C **********************************************************************INT00020
- C THIS IS VERSION 1.57 WHICH USES LESS MEMORY AND IS FASTER THAN *INT00030
- C VERSION 1.55 AND CALC THE SHORT-RANGED DIR CORR FUNCTION(DCS) INT00040
- C PROGRAM RESTRICTED TO SYMMETRICALLY CHARGED ELECTROLYTES. *INT00050
- C SOLVES HNC EQ BY ITERATIONUSING ISEINGS FAST FOURIER TRANSFORM FFT3S*INT00060
- C IF (CONC.EQ.0) READS A NEW SET OF ION AND WELL SIZE PARAMETERS. *INT00070
- C IF REM=0.0, NO REMAINDER ADDED TO G(I,J), T(2,4), G4(I,J) *INT00080
- C IF REM.NE.0.0, DEBYE-HUCKEL REMAINDER ADDED *INT00090
- C PRIMITIVE MODEL PLUS TWO SQUARE WELLS *INT00100
- C SPAC IS NEAREST FRACTION OF A(1,2) THAT IS CLOSEST TO *INT00110
- C MAX(COEFF1*DEBYE/FNN,COEFF2*A(1,2)/FNN) *INT00120
- C DELA = DLN(ACTCO)/DLNC, CALC FROM COMPRESSIBILITY EQUATION. *INT00130
- C DELB = LIMIT(DELB) AS BETA GOES TO ZERO *INT00140
- C BETA MEAS DEVIATION FROM ELECTRONEUTRALITY *INT00150
- C OSM IS CALCULATED FROM THE PRESSURE EQUATION *INT00160
- C XENERGY AND DFX/DP FROM INDEP EQNS ASSUMING NO TEMP AND PRESS DEP *INT00170
- C OF THE SHORT RANGE POTENTIAL XENERGY IS IN CALS PER MOLE OF SOLUTE *INT00180
- C XVOLUME IN CC PER MOLE OF SOLUTE *INT00190
- C DRF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE HNC EQN *INT00200
- C DRDF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE DHLL *INT00210
- C CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12-C(1)+C(2) *INT00220
- C CONVGB = 100*(G(1,2)-G(2,1))/GM12 *INT00230
- C DCS(=ZETA) IS THE SHORT-RANGED DIRECT CORRELATION FUNCTION *INT00240
- C CONVERGENCE CRITERION USEFUL ONLY FOR UNSYMMETRICAL SYSTEMS *INT00250
- C IF EX.EQ.2.0, PRINTS XR,DCS,H *INT00260
- C IF EX.EQ.3.0, PRINTS XR,DCS,H AND CHARGE DENSITY *INT00270
- C IF EX.EQ.0.0, QUITS *INT00280
- C IF CHARGE = 0.0, REPEATS CALC WITH SAME SPACING AS CHARGED SYSTEM *INT00290
- C BUT WITH Z1 = Z2 = 0.0.FICTITIOUS CHARGE DENSITY ALSO CALC *INT00300
- C IF ICC.GT.0, ITERATION BEGINS WITH TAU OF PRECEDING CONC *INT00310
- C IF ICC.LT.0, READ TAU GENERATED IN A PREVIOUS EXECUTION *INT00320
- C IF ICC.EQ.0, TAU CAL IN LAMBDA APPROXIMATION *INT00330
- C WRITES TAU ON NDSK2 IN ALL CASES *INT00340
- C SIX SUBRUTNIES ADJ,TRAP,SIMPSON,FFT3S,SETFT,RSA
- C NIONS is number positive and negative ions 2 for NaCl and ZnSO4
- C **************************DIMENSIONS**********************************INT00350
- INTEGER MNM
- PARAMETER (MNM = 2048)
- REAL*8 DATE,TIME INT00360
- DIMENSION A(2,2)
- DIMENSION A1(2,2)
- DIMENSION B(80)
- DIMENSION BA(2,2)
- DIMENSION B1(2,2)
- DIMENSION C(2)
- DIMENSION CB(2,2)
- DIMENSION CONV(MNM,2,2)
- DIMENSION C1(2,2)
- DIMENSION DC(MNM,2,2)
- DIMENSION DRF(MNM)
- DIMENSION DRDF(MNM)
- DIMENSION E(2,2)
- DIMENSION FNU(2)
- DIMENSION F2(2)
- DIMENSION F4(2)
- DIMENSION G(2,2)
- DIMENSION GA(2,2)
- DIMENSION GB(2,2)
- DIMENSION GC(2,2)
- DIMENSION GM(2,2)
- DIMENSION GW(2,2)
- DIMENSION G4(2,2)
- C DIMENSION H(MNM,2,2)
- REAL H(MNM,2,2)
- DIMENSION LINT(2,2)
- DIMENSION MINT(2,2)
- DIMENSION NINT(2,2)
- DIMENSION PHI(MNM,2,2)
- DIMENSION P1(MNM,2,2)
- DIMENSION Q(MNM,2,2)
- DIMENSION Q1(MNM,2,2)
- DIMENSION RHO(2,MNM)
- DIMENSION RQO(2,MNM)
- DIMENSION R2(2)
- DIMENSION R4(2)
- DIMENSION T(2,2)
- DIMENSION TAU(MNM,2,2)
- DIMENSION TCONV(MNM,2,2)
- DIMENSION TH(MNM,2,2)
- DIMENSION TM(2,2)
- DIMENSION TPHI(MNM,2,2)
- DIMENSION TQ(MNM,2,2)
- DIMENSION TTAU(MNM,2,2)
- DIMENSION TW(2,2)
- DIMENSION TZETA(MNM,2,2)
- DIMENSION W(2,2)
- C DIMENSION XR(MNM)
- REAL XR(MNM)
- DIMENSION Y(2)
- DIMENSION Y1(MNM)
- DIMENSION Z(2)
- DIMENSION ZA(2)
- DIMENSION ZETA(MNM,2,2)
- EQUIVALENCE (H,TH,PHI,TPHI,TTAU)
- EQUIVALENCE (ZETA,TZETA,RHO,DRF,DC)
- EQUIVALENCE (TAU,P1,CONV,TCONV)
- EQUIVALENCE (RQO,DRDF)
- EQUIVALENCE (XR,Y1)
- EQUIVALENCE (B(6),Z(1)),(B(8),D),(B(9),TEMP),(B(10),A1(1,1)),
- *(B(14),COEFF1),(B(15),DEBYE),(B(16),SPAC),(B(17),FNN),(B(18),CONC)
- *,(B(19),DELA),(B(20),OSM),(B(21),CONTAC),(B(22),REMAIN),
- *(B(23),FITER),(B(24),G(1,1)),(B(28),T(1,1)),(B(32),FNGARB),
- *(B(33),DELB),(B(34),DIFF),(B(35),CONVG1),(B(36),CONVG2),
- *(B(37),GA(1,1)),(B(5),DT),(B(4),DP),(B(41),GB(1,1)),
- *(B(45),G4(1,1)),(B(49),F2(1)),(B(51),F4(1)),(B(52),COEFF2)
- COMMON /JAY/ H,XR
- C******************************* Hunting for issues********************
- C Gfortran has some other nifty flags for ID of various issues:
- C -ffpe-trap=invalid,zero,overflow
- C gfortran -g -ffpe-trap=invalid,zero,overflow -std=legacy -fcheck=all
- C******************************* DATE TIME FIX ************************
- CHARACTER*25 date_string
- CHARACTER*100 command
- command = 'date "+%Y-%m-%d %H:%M:%S" > date.txt'
- CALL SYSTEM(command)
- OPEN(UNIT=1, FILE='date.txt', STATUS='OLD')
- READ(1, '(A)') date_string
- CLOSE(1)
- C ******************************* CONSTANTS ****************************INT00600
- C PROGRAM CONSTANTS *INT00610
- DATA PI/3.141593/ INT00620
- PRINT *, date_string
- C CALL TIMER(DATE,TIME,IVIRT,ITOT) INT00630
- NRDR = 5 INT00640
- NPRT = 7 INT00650
- NPUN = 7 INT00660
- NDSK = 8 INT00670
- NDSK2 = 9 INT00680
- NIONS = 2 INT00690
- PI2 = 2.*PI INT00700
- PI4 = 4.*PI INT00710
- PI43 = PI4/3. INT00720
- PIPI = PI*PI INT00730
- C **********************************************************************INT00740
- C CONSTANTS FOR SUBROUTINE FFT3S *INT01430
- F0 = 0. INT01440
- RL = SPAC1*FNN INT01450
- C **********************************************************************INT01460
- C*****used to stop crashing of H MINT PHI TPHI TCONV CONV ZETA*******
- I = 1
- J = 1
- IF (J .LT. 1) THEN
- J = J + 1
- END IF
- IF (I .LT. 1) THEN
- I = I + 1
- END IF
- DO 10 I = 1,80 INT00750
- 10 B(I) = 0.0 INT00760
- B(3) = 1.57 INT00770
- C20 READ(NRDR,9001) M1 INT00780
- C20 READ(NRDR,*) M1 INT00780
- C READ(NRDR,*) (Z(I),FNU(I),I=1,NIONS)
- M1 = 11
- C Z(1) = 1.0
- Z(1) = 2.0
- Z(2) = -1.0 INT00790
- ZE = Z(1) INT00800
- ZA(1) = Z(1) INT00810
- ZA(2) = Z(2) INT00820
- FNU(1) = 1.0
- FNU(2) = 2.0
- C25 READ(NRDR,*) D,TEMP,DT,DP
- D = 78.358
- T = 298.15
- TEMP = 298.15
- DT = 0.0
- DP = 0.0
- C30 READ(NRDR,*) A(1,1),BA(1,1),CB(1,1),A(1,2),BA(1,2),CB(1,2),A(2, INT00840
- C *2),BA(2,2),CB(2,2) INT00850
- A(1,1) = 4.2
- A(1,2) = 4.2
- A(2,1) = A(1,2) INT00860
- A(2,2) = A (1,1)
- BA(1,1) = 0
- BA(2,2) = BA(1,1)
- BA(1,2) = 0.0
- BA(2,1) = BA(1,2)
- CB(1,1) = 0.0
- CB(2,2) = CB(1,1)
- CB(1,2) = 0.0 INT00870
- CB(2,1) = CB(1,2) INT00880
- C35 READ(NRDR,*) E(1,1),W(1,1),E(1,2),W(1,2),E(2,2),W(2,2) INT00890
- E(1,1)= 0.0
- E(1,2)= 0.0
- E(2,1) = E(1,2) INT00900
- E(2,2) = E(1,1) INT00900
- W(1,1) = 0.0
- W(1,2) = 0.0
- W(2,1) = W(1,2)
- W(2,2) = W(1,1)
- C40 READ(NRDR,*)CONC,COEFF1,COEFF2,EX,REM,CHARGE,MAX,I2,NCD,IC, INT00920
- C *NR,IR,ICC
- CONC = 0.1 INT00930
- COEFF1=1.0
- COEFF2=1.0
- EX = 3.0
- REM = 1.0
- CHARGE = 1.0
- MAX = 6
- I2 =2
- ICC= 1
- C NCD =
- C IC =
- IR = 1
- NR = 5
- C REWIND NDSK2 INT00940
- C IF(ICC.LT.0)READ(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=INT00950
- C *1,N) INT00960
- C REWIND NDSK2 INT00970
- IF (EX.EQ.0.0) GO TO 5050 INT00980
- C IF(CONC.EQ.0.0) GO TO 30 INT00990
- 9001 FORMAT(I3) INT01000
- 9002 FORMAT(8F10.5) INT01010
- 9003 FORMAT(4F10.5) INT01020
- 9004 FORMAT(9F5.2) INT01030
- 9005 FORMAT(6F10.5) INT01040
- 9006 FORMAT(1F10.7,5F5.2,7I5) INT01050
- WRITE(NPRT,9100)B(3), date_string
- 9100 FORMAT(1H1, 38X, 'HNC VERSION ', F4.2, ' - IBM 370/148 - ', A25)
- C WRITE(NPRT,9100)B(3),DATE,TIME INT01060
- C9100 FORMAT(1H1,38X,'HNC VERSION ',F4.2,' - IBM 370/148 - ',A8,2X,A8) INT01070
- C**** N=2**M1=2048 for M1=11, ***************************************
- N=2**M1 INT01090
- N2=N/2 INT01100
- FNN = N INT01120
- FC = 0.0 INT01130
- FNO = 0.0 INT01140
- CC = 0.0 INT01150
- DO 50 I = 1,NIONS INT01160
- FNO = FNU(I)+FNO INT01170
- C(I) = FNU(I)*CONC*6.0233*1.E-04 INT01180
- CC = C(I)+CC INT01190
- 50 FC = C(I)*Z(I)**2+FC
- C********************************************************************** INT01200
- C Electron charge e=4.8032x10E10 cm**3/2 g**1/2 s**-2 (esu)
- C Boltzconst k=1.3803x10E-16 cm**2 s**2/K INT01200
- C Bjerrum length |z+z+|e**2/DkT=16.714x1.E+04/DT
- C***********************************************************************
- FL1 = 16.7094*1.E+04/(D*TEMP) INT01210
- FL = PI4*FL1 INT01220
- FK = (FC*FL)**0.5 INT01230
- FKFK = FK*FK INT01240
- DEBYE = 1./FK INT01250
- FKA = A(1,2)/DEBYE INT01260
- C REWIND NDSK INT01270
- C WRITE(NDSK,9101)A(1,1),A(1,2),A(2,2),D,TEMP,CONC,FKA,DEBYE,Z(1),Z(INT01280
- C *2) INT01290
- 9101 FORMAT(3F8.5,2F8.4,3F10.6,2F5.1) INT01300
- SPAC1 = COEFF1*DEBYE/FNN INT01310
- SPAC2 = COEFF2*A(1,2)/FNN INT01320
- IF(SPAC2.GT.SPAC1) SPAC1=SPAC2 INT01330
- S1 = A(1,2)/SPAC1 INT01340
- L = S1+1. INT01350
- S2 = L INT01360
- IF (S2-S1.GT.0.5) GO TO 60 INT01370
- SPAC1 = A(1,2)/S2 INT01380
- GO TO 70 INT01390
- 60 SPAC1 = A(1,2)/(S2-1.) INT01400
- 70 CONTINUE INT01410
- C **********************************************************************INT01420
- C CONSTANTS FOR SUBROUTINE FFT3S *INT01430
- F0 = 0. INT01440
- RL = SPAC1*FNN INT01450
- C FX,FY and FZ are Debye Huckel remainders for ?
- C **********************************************************************INT01460
- FI = 1.0 INT01470
- SPAC = SPAC1 INT01480
- DO 80 K = 1,N INT01490
- XR(K) = SPAC1*FI INT01500
- 80 FI = FI+1. INT01510
- XN = XR(N) INT01520
- FKN = FK*XN INT01530
- IF (REM.EQ.0.0) GO TO 90 INT01540
- FX = FL*(1.+FKN)*EXP(-FKN)/FKFK INT01550
- FY = FL*EXP(-FKN)/FK INT01560
- FZ = FL*EXP(-FKN)*(FKN**3+3.*FKN**2+6.*FKN+6.)/(FKFK*FKFK) INT01570
- GO TO 100 INT01580
- 90 FX = 0.0 INT01590
- FY = 0.0 INT01600
- FZ = 0.0 INT01610
- 100 CONTINUE INT01620
- DO 110 I = 1,NIONS INT01630
- DO 110 J = 1,NIONS INT01640
- CALL ADJ (A(I,J),A1(I,J),MINT(I,J),SPAC1) INT01650
- CALL ADJ (BA(I,J),B1(I,J),NINT(I,J),SPAC1) INT01660
- CALL ADJ (CB(I,J),C1(I,J),LINT(I,J),SPAC1) INT01670
- 110 CONTINUE INT01680
- GO TO 130 INT01690
- 120 PRINT *, date_string
- C 120 CALL TIMER(DATE,TIME,IVIRT,ITOT) INT01700
- WRITE(NPRT,9100)B(3),date_string INT01710
- C WRITE(NPRT,9100)B(3),DATE,TIME INT01710
- Z(1) = 0.0 INT01720
- Z(2) = 0.0 INT01730
- ZE = 0.0 INT01740
- CHARGE = 20.0 INT01750
- REM = 0.0 INT01760
- DEBYE = 0.0 INT01770
- BETA = 0.0 INT01780
- CONVGA = 0.0 INT01790
- 130 WRITE(NPRT,9103) Z(1),Z(2),D,TEMP,DT,DP,A1(1,1),A1(2,2),A1(1,2), INT01800
- *B1(1,1),B1(2,2),B1(1,2),C1(1,1),C1(2,2),C1(1,2),E(1,1),E(2,2), INT01810
- *E(1,2),W(1,1),W(2,2),W(1,2),COEFF1,COEFF2,SPAC1,N,REM,CONC,DEBYE INT01820
- 9103 FORMAT('Z1=',F5.1,3X,'Z2=',F5.1,/,'DIEL CONST=',F7.3,3X,'TEMP=', INT01830
- *F7.3,'DLND/DLNT=',F10.5,3X,'DLND/DP=',F7.3,/,'A11=',F6.3,3X,'A22='INT01840
- *,F6.3,3X,'A12=',F6.3,/,'B11=',F6.3,3X,'B22=',F6.3,3X,'B12=',F6.3,/INT01850
- *'C11=',F6.3,3X,'C22=',F6.3,3X,'C12=',F6.3,/, INT01860
- *'E11=',F6.3,3X,'E22=',F6.3,3X,'E12=',F6.3,/, INT01870
- *'W11=',F6.3,3X,'W22=',F6.3,3X,'W12=',F6.3,/, INT01880
- *'COEFF1=',F7.4,3X,'COEFF2=',F7.4,3X,'SPAC1=',F7.4,3X,'N=',I4,3X, INT01890
- *'REM=',F7.4,/,'CONC=',F10.6,3X,'DEBYE=',F10.5,/) INT01900
- C *** CALC OF Q,Q1,AND TQ INT01910
- C WRITE(6,5558) INT01920
- 5558 FORMAT(1X,'HO') INT01930
- DO 180 I = 1,NIONS INT01940
- DO 180 J = 1,NIONS INT01950
- AB = A1(I,J)+B1(I,J) INT01960
- BC = AB+C1(I,J) INT01970
- DO 180 K = 1,N INT01980
- Q(K,I,J) = -Z(I)*Z(J)*FL*EXP(-FK*XR(K))/(PI4*XR(K)) INT01990
- IF (XR(K) .LT. A1(I,J)) GO TO 170 INT02000
- IF(ABS(XR(K)-AB).LT.1.E-04) GO TO 160 INT02010
- IF ((XR(K)-AB).LT.1.E-07) GO TO 150 INT02020
- IF(ABS(XR(K)-BC).LT.1.E-04) GO TO 140 INT02030
- IF((XR(K)-BC).LT.1.E-07) GO TO 160 INT02040
- 140 Q1(K,I,J) = EXP(Q(K,I,J)) INT02050
- GO TO 180 INT02060
- 150 Q1(K,I,J) = EXP(Q(K,I,J)-E(I,J)) INT02070
- GO TO 180 INT02080
- 160 Q1(K,I,J) = EXP(Q(K,I,J)-W(I,J)) INT02090
- GO TO 180 INT02100
- 170 Q1(K,I,J) = 0.0 INT02110
- 180 CONTINUE INT02120
- C WRITE(6,5557) INT02130
- 5557 FORMAT(1X,'HA') INT02140
- DO 190 I = 1,NIONS INT02150
- DO 190 J = 1,NIONS INT02160
- 190 CALL FFT3S(F0,Q(1,I,J),M1,RL,TQ(1,I,J),1) INT02170
- C***************** Distributution ************************************
- C CALC OF APPROX DISTRIBUTION FUNCTIONS after counter I1 =0 *INT02190
- 1010 I1 = 0 INT02210
- IF (ICC.NE.0) GO TO 1070 INT02220
- DO 1020 I = 1,NIONS INT02230
- DO 1020 J = 1,NIONS INT02240
- DO 1020 K = 1,N INT02250
- 1020 PHI(K,I,J) = Q1(K,I,J)-1. INT02260
- DO 1030 I = 1,NIONS INT02270
- DO 1030 J = 1,NIONS INT02280
- PHI(MINT(I,J),I,J) = (PHI(MINT(I,J),I,J)-1.)*.5 INT02290
- 1030 CALL FFT3S(F0,PHI(1,I,J),M1,RL,TPHI(1,I,J),1) INT02300
- DO 1040 I = 1,NIONS INT02310
- DO 1040 J = 1,NIONS INT02320
- DO 1040 K = 1,N INT02330
- TCONV(K,I,J) = 0.0 INT02340
- DO 1040 L = 1,NIONS INT02350
- 1040 TCONV(K,I,J) = TPHI(K,I,L)*TPHI(K,L,J)*C(L)+TCONV(K,I,J) INT02360
- DO 1050 I = 1,NIONS INT02370
- DO 1050 J = 1,NIONS INT02380
- 1050 CALL FFT3S(F0,TCONV(1,I,J),M1,RL,CONV(1,I,J),-1) INT02390
- DO 1060 I = 1,NIONS INT02400
- DO 1060 J = 1,NIONS INT02410
- DO 1060 K = 1,N INT02420
- 1060 TAU(K,I,J) = CONV(K,I,J)+Q(K,I,J)*XR(K)*FK*.5 INT02430
- 1070 DO 1080 I = 1,NIONS INT02440
- DO 1080 J = 1,NIONS INT02450
- DO 1075 K = 1,N INT02460
- 1075 H(K,I,J) = Q1(K,I,J)*(1.+TAU(K,I,J))-1. INT02470
- IF (MINT(I,J) < 1) THEN
- MINT(I,J) = 1
- END IF
- H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT02480
- DO 1080 K = 1,N INT02490
- 1080 ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J) INT02500
- IF (ICC.NE.0) GO TO 1090 INT02510
- WRITE(NPRT,9104) INT02520
- 9104 FORMAT(38H0RESULTS OF TRIANGLE TERM IN DIST FUNC) INT02530
- GO TO 1100 INT02540
- 1090 IF(ICC.GT.0) GO TO 1112 INT02550
- WRITE(NPRT,9105) INT02560
- 9105 FORMAT(1H0,'RESULTS USING PRECEDING TAU AT PRESENT CONC') INT02570
- GO TO 1100 INT02580
- 1112 WRITE(NPRT,9102) INT02590
- 9102 FORMAT(1H0,'RESULTS USING PREVIOUS TAU STORED ON DISK') INT02600
- 1100 WRITE(NPRT,9106) INT02610
- 9106 FORMAT(6X,4HDELA,9X,4HDELB,9X,4HDIFF,10X,3HOSM,7X,7HXENERGY,5X, INT02620
- *7HXVOLUME,4X,4HITER,6X,4HBETA,7X,6HCONVGA,7X,6HCONVGB) INT02630
- GO TO 3010 INT02640
- C***************** HNC ***********************************************
- C SOLVE HNC INTEGRAL EQUATION BY ITERATION *INT02660
- 2010 WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG INT02680
- WRITE(NPRT,5555)ZE INT02690
- IF(ZE.EQ.0.0) GO TO 2020 INT02700
- 5555 FORMAT(1X,F10.5) INT02710
- WRITE(NPRT,9108) INT02720
- WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS) INT02730
- WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4 INT02740
- 2020 WRITE(NPRT,9111) INT02750
- 9111 FORMAT(19H0RESULTS OF HNC EQN) INT02760
- WRITE(NPRT,9106) INT02770
- 2030 CONTINUE INT02780
- DO 2040 I = 1,NIONS INT02790
- DO 2040 J = 1,NIONS INT02800
- 2040 CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1) INT02810
- DO 2050 I = 1,NIONS INT02820
- DO 2050 J = 1,NIONS INT02830
- 2050 CALL FFT3S(F0,ZETA(1,I,J),M1,RL,TZETA(1,I,J),1) INT02840
- DO 2070 I = 1,NIONS INT02850
- DO 2070 J = 1,NIONS INT02860
- DO 2070 K = 1,N INT02870
- SUM1 = 0.0 INT02880
- DO 2060 L = 1,NIONS INT02890
- 2060 SUM1 = TZETA(K,I,L)*TH(K,L,J)*C(L)+SUM1 INT02900
- 2070 P1(K,I,J) = SUM1 INT02910
- DO 2090 I = 1,NIONS INT02920
- DO 2090 J = 1,NIONS INT02930
- DO 2090 K = 1,N INT02940
- P2 = 0.0 INT02950
- P3 = 0.0 INT02960
- DO 2080 L = 1,NIONS INT02970
- P2 = TQ(K,I,L)*TZETA(K,L,J)*C(L)+P2 INT02980
- 2080 P3 = TQ(K,I,L)*P1(K,L,J)*C(L)+P3 INT02990
- 2090 TTAU(K,I,J) = P1(K,I,J)+P2+P3 INT03000
- DO 2100 I = 1,NIONS INT03010
- DO 2100 J = 1,NIONS INT03020
- 2100 CALL FFT3S(F0,TTAU(1,I,J),M1,RL,TAU(1,I,J),-1) INT03030
- I1 = I1+1 INT03040
- 2110 DO 2120 I = 1,NIONS INT03050
- DO 2120 J = 1,NIONS INT03060
- C WRITE(7,5556) INT03070
- 5556 FORMAT(1X,'HI') INT03080
- DO 2115 K = 1,N INT03090
- C somethign about this raises the floating point exception
- C looking deeper, the math library doesnt like that something here has
- C an inf value. if is a valid float though. line 39 of
- C math_errf.c library.
- C it probably would preffer IMPLICIT NONE
- C
- 2115 H(K,I,J) = Q1(K,I,J)*EXP(TAU(K,I,J))-1. INT03100
- H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT03110
- DO 2120 K = 1,N INT03120
- 2120 ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J) INT03130
- IF (I1.EQ.MAX) GO TO 3010 INT03140
- IF ((I1/I2)*I2.EQ.I1) GO TO 3010 INT03150
- GO TO 2030
- C**************** XEnergy*********************************************
- C DELA,DELB,OSM,XENERGY,DFX/DP *INT03180
- 3010 CONTAC = 0.0 INT03200
- SQ1 = 0.0 INT03210
- SQ2 = 0.0 INT03220
- REMAIN = 0.0 INT03230
- XE = 0.0 INT03240
- DO 3020 I = 1,NIONS INT03250
- DO 3020 J = 1,NIONS INT03260
- NA = MINT(I,J) INT03270
- HNA = H(NA,I,J) INT03280
- H(NA,I,J) = HNA*2.+1. INT03290
- AA = A(I,J)**3 INT03300
- NB = NA +NINT(I,J) INT03310
- BB = (A(I,J)+BA(I,J))**3 INT03320
- NC = NB+LINT(I,J) INT03330
- AC = (A(I,J)+BA(I,J)+CB(I,J))**3 INT03340
- AD = PI4*(A1(I,J)-A(I,J)) INT03350
- CALL TRAP(NA,NB,I,J,SUM1,SUM2,SUM3) INT03360
- T(I,J) = PI4*SUM1*SPAC1 INT03370
- G(I,J) = PI4*SUM2*SPAC1 INT03380
- G4(I,J) = PI4*SUM3*SPAC1 INT03390
- GM(I,J) = G(I,J) INT03400
- TM(I,J) = T(I,J) INT03410
- XE = XE + FNU(I)*FNU(J)*E(I,J)*(PI43*(BB-AA)+G(I,J)) INT03420
- CALL TRAP(NB,NC,I,J,SUM1,SUM2,SUM3) INT03430
- T(I,J) = PI4*SUM1*SPAC1 INT03440
- G(I,J) = PI4*SUM2*SPAC1 INT03450
- G4(I,J) = PI4*SUM3*SPAC1+G4(I,J) INT03460
- GW(I,J) = G(I,J) INT03470
- TW(I,J) = T(I,J) INT03480
- XE = XE+FNU(I)*FNU(J)*W(I,J)*(PI43*(AC-BB)+G(I,J)) INT03490
- ERROR1 = -Z(I)*Z(J)*FY+AD*H(NA,I,J)*XR(NA) INT03500
- ERROR2 = -Z(I)*Z(J)*FX+AD*H(NA,I,J)*XR(NA)**2 INT03510
- ERROR3 = -Z(I)*Z(J)*FZ+AD*H(NA,I,J)*XR(NA)**4 INT03520
- C PRINT *, "*****From xEnergy*******"
- C PRINT *, "NC",NC,"N",N,"I",I,"J",J
- C PRINT *,"SUM1",SUM1,"SUM2",SUM2,"SUM3",SUM3
- CALL SIMP(NC,N,I,J,SUM1,SUM2,SUM3) INT03530
- T(I,J) = PI43*SPAC1*SUM1-PI2*A(I,J)**2+TW(I,J)+TM(I,J)+ERROR1 INT03540
- G(I,J) = PI43*(SPAC1*SUM2-A(I,J)**3)+GM(I,J)+GW(I,J)+ERROR2 INT03560
- G4(I,J) = PI43*SPAC1*SUM3-2.513234*A(I,J)**5+G4(I,J)+ERROR3 INT03570
- GA(I,J) = H(NA,I,J)+1. INT03580
- GB(I,J) = H(NB,I,J)+1. INT03590
- GC(I,J) = H(NC,I,J)+1. INT03600
- CONTAC = CONTAC + FNU(I)*FNU(J)*GA(I,J)*AA INT03610
- SQ1 = SQ1+FNU(I)*FNU(J)*(GB(I,J)*BB*(1.-EXP(-E(I,J)+W(I,J)))) INT03620
- SQ2 = SQ2+FNU(I)*FNU(J)*(GC(I,J)*AC*(1.-EXP(-W(I,J)))) INT03630
- REMAIN = REMAIN + Z(I)*Z(J)*FNU(I)*FNU(J)*T(I,J) INT03640
- IF(I1.EQ.MAX) GO TO 3011 INT03650
- H(NA,I,J) = HNA INT03660
- GO TO 3020 INT03670
- 3011 DC(NA,I,J) = DC(NA,I,J)+HNA+1. INT03680
- 3020 CONTINUE INT03690
- CONTAC = CONTAC*2.09439*CC/FNO**2 INT03700
- SQ1 = SQ1*2.09439*CC/FNO**2 INT03710
- SQ2 = SQ2*2.09439*CC/FNO**2 INT03720
- REMAIN = REMAIN*FL1*CC/(6.*FNO**2) INT03730
- OSM = 1.+CONTAC+SQ1+SQ2+REMAIN INT03740
- XE = 1.97827*TEMP*(3.*(1.+DT)*REMAIN*FNO+XE*CC*0.5/FNO) INT03750
- XFP = -8.31467*1.E-05*TEMP*DP*3.*REMAIN*FNO INT03760
- GM12 = (G(1,2)+G(2,1))*.5 INT03770
- CONVGB = 100*(G(1,2)-G(2,1))/GM12 INT03780
- IF (ZE.EQ.0.0) GO TO 3050 INT03790
- S2 = 0.0 INT03800
- S4 = 0.0 INT03810
- DO 3040 I = 1,NIONS INT03820
- SUM1 = 0.0 INT03830
- SUM2 = 0.0 INT03840
- DO 3030 J = 1,NIONS INT03850
- SUM1 = SUM1 + C(J)*Z(J)*G(I,J) INT03860
- 3030 SUM2 = SUM2 + C(J)*Z(J)*G4(I,J) INT03870
- F2(I) = SUM1 INT03880
- F4(I) = SUM2 INT03890
- Y(I) = 6.*Z(I)/FKFK INT03900
- R2(I) = 1.+F2(I)/Z(I) INT03910
- R4(I) = 1. + F4(I)/Y(I) INT03920
- S2 = C(I)*Z(I)*F2(I)+S2 INT03930
- 3040 S4 = C(I)*Z(I)*F4(I)+S4 INT03940
- CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12 - C(1)+C(2) INT03950
- CS2 = FKFK/FL INT03960
- CS4 = 6./FL INT03970
- E2 = 1.+S2/CS2 INT03980
- E4 = 1.+S4/CS4 INT03990
- BETA = FNO*R2(1)*R2(2)/(FNU(1)*R2(2)+FNU(2)*R2(1)) INT04000
- 3050 CONTINUE INT04010
- FG = 1.0/(CC*GM12) INT04020
- DELA = (FG-1.0-FG*BETA)/(1.0+FG*BETA) INT04030
- DELB = FG-1.0 INT04040
- DIFF = DELA-DELB INT04050
- ITER = I1 INT04060
- FITER = ITER INT04070
- WRITE(NPRT,9112) DELA,DELB,DIFF,OSM,XE,XFP,ITER,BETA,CONVGA,CONVGBINT04080
- 9112 FORMAT( 6E13.4,I5,3E13.4) INT04090
- IF (I1.EQ.MAX) GO TO 3060 INT04100
- IF(I1.EQ.0) GO TO 2010 INT04110
- IF ((I1/I2)*I2.EQ.I1) GO TO 2030 INT04120
- 3060 WRITE(NPRT,9113) INT04130
- 9113 FORMAT(15H0MATRIX G(I,J),21X,15HMATRIX GM(I,J),20X, INT04140
- *15HMATRIX GW(I,J),20X,15HMATRIX G4(I,J)) INT04150
- WRITE(NPRT,9114) ((G(I,J),J=1,NIONS),(GM(I,L),L=1,NIONS),(GW(I,M),INT04160
- *M=1,NIONS),(G4(I,K),K=1,NIONS),I=1,NIONS) INT04170
- 9114 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6,6X,2E14.6) INT04180
- WRITE(NPRT,9115) INT04190
- 9115 FORMAT(15H0MATRIX T(I,J),21X,15HMATRIX TM(I,J),20X, INT04200
- *15HMATRIX TW(I,J)) INT04210
- WRITE(NPRT,9116)((T(I,J),J=1,NIONS),(TM(I,L),L=1,NIONS),(TW(I,M),MINT04220
- *=1,NIONS),I=1,NIONS) INT04230
- 9116 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6) INT04240
- WRITE(NPRT,9117) INT04250
- 9117 FORMAT(16H0MATRIX GA(I,J),20X,15HMATRIX GB(I,J),20X, INT04260
- *15HMATRIX GC(I,J)) INT04270
- WRITE(NPRT,9118) ((GA(I,J),J=1,NIONS),(GB(I,L),L=1,NIONS),(GC(I,M)INT04280
- *,M=1,NIONS),I=1,NIONS) INT04290
- 9118 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6) INT04300
- WRITE(NPRT,9119) INT04310
- 9119 FORMAT(17H0MATRIX MINT(I,J),3X,16HMATRIX NINT(I,J),3X, INT04320
- *16HMATRIX LINT(I,J)) INT04330
- WRITE(NPRT,9120) ((MINT(I,J),J=1,NIONS),(NINT(I,L),L=1,NIONS),(LININT04340
- *T(I,M),M=1,NIONS),I=1,NIONS) INT04350
- 9120 FORMAT(1X,2I4,13X,2I4,13X,2I4) INT04360
- WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG INT04370
- 9107 FORMAT(8H0CONTAC=,1P,E12.4,3X,'SQUARE1=',1P,E12.4,3X,'SQUARE2=', INT04380
- *1P,E12.4,3X,'REMAIN=',1P,E12.4,3X,'GAMMA=',1P,E12.4) INT04390
- IF (ZE.EQ.0.0) GO TO 3070 INT04400
- WRITE(NPRT,9108) INT04410
- 9108 FORMAT(6X,1HZ,13X,2HF2,12X,2HR2,12X,1HY,13X,2HF4,12X,2HR4) INT04420
- WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS) INT04430
- 9109 FORMAT(1X,6E14.6) INT04440
- WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4 INT04450
- 9110 FORMAT('S2=',E14.6,2X,'CS2=',1P,E14.6,2X,'E2=',1P,E14.6,2X,'S4=' INT04460
- *,1P,E14.6,2X,'CS4=',1P,E14.6,2X,'E4=',1P,E14.6) INT04470
- 3070 CONTINUE INT04480
- WRITE(NPRT,9121) INT04490
- 9121 FORMAT(7X,5HXR(N),7X,9H Q(1,1,N),5X,8HH(1,1,N),5X,9H Q(1,2,N), INT04500
- *6X,8HH(1,2,N),5X,9H Q(2,2,N),6X,8HH(2,2,N)) INT04510
- WRITE(NPRT,9122) XR(N),Q(N,1,1),H(N,1,1),Q(N,1,2),H(N,1,2), INT04520
- *Q(N,2,2),H(N,2,2) INT04530
- 9122 FORMAT(7E14.4//) INT04540
- IF(EX.NE.2.0.AND.EX.NE.3.0) GO TO 4070 INT04550
- 3080 WRITE(NPRT,9123) INT04560
- 9123 FORMAT(7X,5HXR(K),6X,10H DC(1,1,K),5X,8HH(1,1,K),5X,10H DC(1,2,K),INT04570
- *5X,8HH(1,2,K),5X,10H DC(2,1,K),5X,8HH(2,1,K),5X,10H DC(2,2,K), INT04580
- *5X,8HH(2,2,K)) INT04590
- WRITE(NPRT,9124)(XR(K),((DC(K,I,J),H(K,I,J),J=1,NIONS),I=1,NIONS) INT04600
- *,K=1,N) INT04610
- WRITE(NDSK,9125)C1,C2 INT04620
- WRITE(NDSK,9133)FO,N,M1,RL,((MINT(I,J),I=1,2),J=1,2) INT04630
- 9133 FORMAT(F8.4,2I4,F8.4,4I15.6) INT04640
- WRITE(NDSK,9126)(XR(K),H(K,1,1),H(K,1,2),K=1,N) INT04650
- WRITE(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=1,N) INT04660
- REWIND NDSK2 INT04670
- 9124 FORMAT(1X,1P,9E14.4) INT04680
- C ERIC some issue with this line format
- 9125 FORMAT(1X,A8,2E16.6) INT04690
- 9126 FORMAT(2E16.7) INT04700
- IF (ZE.EQ.0.0) GO TO 4030 INT04710
- IF (EX.NE.3.0) GO TO 4070 INT04720
- C **********************************************************************INT04730
- C CALC OF CHARGE DENSITY *INT04740
- C **********************************************************************INT04750
- DO 4020 K = 1,N INT04760
- DO 4020 I = 1,NIONS INT04770
- SUM1 = 0.0 INT04780
- SUM2 = 0.0 INT04790
- DO 4010 J = 1,NIONS INT04800
- SUM1 = SUM1 + FNU(J)*Z(J)*H(K,I,J) INT04810
- 4010 SUM2 = SUM2 + FNU(J)*Z(J)*Q(K,I,J) INT04820
- RHO(I,K) = SUM1*XR(K)**2 INT04830
- 4020 RQO(I,K) = SUM2*XR(K)**2 INT04840
- GO TO 4060 INT04850
- 4030 Z(1) = 1.0 INT04860
- Z(2) = -1.0 INT04870
- DO 4050 K = 1,N INT04880
- DO 4050 I = 1,NIONS INT04890
- SUM1 = 0.0 INT04900
- DO 4040 J = 1,NIONS INT04910
- 4040 SUM1 = SUM1+FNU(J)*Z(J)*H(K,I,J) INT04920
- RHO(I,K) = SUM1*XR(K)**2 INT04930
- 4050 RQO(I,K) = 0.0 INT04940
- 4060 WRITE(NPRT,9127) INT04950
- 9127 FORMAT('1CHARGE DENSITY'///1X,9X,'XR',11X,'RHO1',10X,'RQO1',11X, INT04960
- *'RHO2',10X,'RQO2') INT04970
- C PRINT *, 'NCD:' , NCD
- C PRINT *, 'K:' , K
- C PRINT *, 'IC:' , IC
- WRITE(NPRT,9128)(XR(K),RHO(1,K),RQO(1,K),RHO(2,K),RQO(2,K),K=1,NCDINT04980
- *,IC) INT04990
- PRINT *, 'NR:' , NR
- 9128 FORMAT(1X,1P,5E14.4) INT05000
- WRITE(NDSK,9126)(RHO(1,K),RHO(2,K),K=1,N2) INT05010
- 9129 FORMAT(2E16.7) INT05020
- REWIND NDSK INT05030
- 4070 CONTINUE INT05040
- B(1) = 0.0 INT05050
- B(2) = 1.0 INT05060
- Z(1) = ZA(1) INT05070
- Z(2) = ZA(2) INT05080
- ZE = Z(1) INT05090
- IF (ZE.EQ.0.0) GO TO 5040 INT05100
- C **********************************************************************INT05110
- C CALC OF DIELECTRIC RESPONSE FUNCTION INT05120
- C **********************************************************************INT05130
- DO 5010 I = 1,NIONS INT05140
- DO 5010 J = 1,NIONS INT05150
- H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT05160
- 5010 CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1) INT05170
- SPACT = PI2/(2.*FNN-1.)/SPAC1 INT05180
- TK = 1.0 INT05190
- DO 5030 K = 1,N INT05200
- Y1(K) = TK*SPACT INT05210
- SUM = 0.0 INT05220
- DO 5020 I = 1,NIONS INT05230
- DO 5020 J = 1,NIONS INT05240
- 5020 SUM = SUM + C(I)*C(J)*Z(I)*Z(J)*TH(K,I,J) INT05250
- Y1K2 = Y1(K)**2 INT05260
- DRF(K) = 1.-(FKFK/Y1K2)*(1.+SUM/FC) INT05270
- DRDF(K) = 1.-FKFK/(FKFK+Y1K2) INT05280
- 5030 TK = TK+1. INT05290
- WRITE(NPRT,9130) INT05300
- 9130 FORMAT('1DIELECTRIC RESPONSE FUNCTION'///1X,3(9X,'YK',12X,'DRF', INT05310
- *11X,'DRDF')) INT05320
- WRITE(NPRT,9131)(Y1(K),DRF(K),DRDF(K),K=1,NR,IR) INT05330
- 9131 FORMAT(1X,1P,9E14.4) INT05340
- 5040 CONTINUE
- C5040 CALL TIMER(DATE,TIME,IVIRT,ITOT) INT05350
- C RVIRT=FLOAT(IVIRT)*1.E-6 INT05360
- C WRITE(NPRT,9132)TIME,DATE,RVIRT INT05370
- WRITE(NPRT,9132)date_string INT05370
- C9132 FORMAT(//////////'END OF EXECUTION AT',A8,' ON ',A8,5X,F8.3, INT05380
- C *' SECONDS USED') INT05390
- 9132 FORMAT(//////////'END EXECUTION ON ',A20) INT05390
- IF (ABS(CHARGE-0.0).LT.1.E-04) GO TO 120 INT05400
- C GOTO 40 INT05410
- C Program loop
- 5050 CONTINUE INT05420
- STOP INT05430
- END INT05440
- C ****Program stops and ends ****************************************
- C************************ SUBROUTINES *******************************
- C************************ ADJ ***************************************
- SUBROUTINE ADJ(A,A1,M,SPAC1) INT05450
- S1 = A/SPAC1 INT05460
- L = S1+1 INT05470
- S2 = L INT05480
- IF(S2-S1.GE.0.5) GO TO 10 INT05490
- A1 = S2*SPAC1 INT05500
- M = L INT05510
- GO TO 20 INT05520
- 10 A1 = (S2-1.)*SPAC1 INT05530
- M = L-1 INT05540
- 20 CONTINUE INT05550
- RETURN INT05560
- END
- C************************ TRAP ***************************************
- SUBROUTINE TRAP(N1,N11,I,J,SUM1,SUM2,SUM3) INT05580
- DIMENSION H(2048,2,2),XR(2048) INT05590
- COMMON/JAY/ H,XR INT05600
- IF (N11.EQ.N1) GO TO 20 INT05610
- SUM1 = (H(N1,I,J)*XR(N1)+H(N11,I,J)*XR(N11))*.5 INT05620
- SUM2 = (H(N1,I,J)*XR(N1)**2+H(N11,I,J)*XR(N11)**2)*.5 INT05630
- SUM3 = (H(N1,I,J)*XR(N1)**4+H(N11,I,J)*XR(N11)**4)*.5 INT05640
- L1 = N1+1 INT05650
- L2 = N11-1 INT05660
- IF (N11.EQ.L1) GO TO 30 INT05670
- DO 10 K = L1,L2 INT05680
- SUM1 = SUM1 + H(K,I,J)*XR(K) INT05690
- SUM3 = SUM3 + H(K,I,J)*XR(K)**4 INT05700
- 10 SUM2 = SUM2 + H(K,I,J)*XR(K)**2 INT05710
- GO TO 30 INT05720
- 20 SUM1 = 0.0 INT05730
- SUM2 = 0.0 INT05740
- SUM3 = 0.0 INT05750
- 30 RETURN INT05760
- END
- C********************** SIMPSON **************************************
- SUBROUTINE SIMP(NC,N,I,J,SUM1,SUM2,SUM3) INT05780
- DIMENSION H(2048,2,2),XR(2048) INT05790
- COMMON /JAY/ H,XR INT05800
- SUM1 = H(NC,I,J)*XR(NC) INT05810
- SUM2 = SUM1*XR(NC) INT05820
- SUM3 = SUM2*XR(NC)**2 INT05830
- C PRINT *, 'NC:',NC
- C Setting NC to 1 for testing??
- NC = 1
- N2 = (N-NC-2)/2
- C PRINT *, '**N2:', N2
- DO 10 K = 1,N2 INT05850
- K4 = NC+2*K-1 INT05860
- K5 = K4+1 INT05870
- K6 = NC+2*N2+1 INT05880
- K7 = K6+1 INT05890
- SUM1 = SUM1 + 4.*H(K4,I,J)*XR(K4)+2.*H(K5,I,J)*XR(K5) INT05900
- SUM3 = SUM3 + 4.*H(K4,I,J)*XR(K4)**4+2.*H(K5,I,J)*XR(K5)**4 INT05910
- 10 SUM2 = SUM2 + 4.*H(K4,I,J)*XR(K4)**2+2.*H(K5,I,J)*XR(K5)**2 INT05920
- SUM1 = SUM1 + 4.*H(K6,I,J)*XR(K6)+2.*H(K7,I,J)*XR(K7) INT05930
- SUM2 = SUM2 + 4.*H(K6,I,J)*XR(K6)**2+2.*H(K7,I,J)*XR(K7)**2 INT05940
- SUM3 = SUM3 + 4.*H(K6,I,J)*XR(K6)**4+2.*H(K7,I,J)*XR(K7)**4 INT05950
- RETURN INT05960
- END
- C**************** Fast Fourier Transform *************************
- SUBROUTINE FFT3S(F0,F,MM,RL,T,MODE) INT05990
- C FAST FOURIER TRANSFORM OF 3-DIMENSIONAL CENTROSYMMETRIC FUNCTION INT06000
- DIMENSION F(2048),T(1) INT06010
- COMMON /FDATA/ N,N2,FF,RTTWO,NC,NS,ND,NR INT06020
- COMMON /FWORK/ W(12288) INT06030
- C DIMENSION OF W SHOULD BE .GE. 3*N INT06040
- C FIXED DATA PREPARATION INT06050
- DATA PI/3.141593/,M/-1/ INT06060
- MD = MM+1 INT06070
- CALL SETFT(MD,M) INT06080
- N = 2**MM INT06090
- RN = FLOAT(N) INT06100
- DR = RL/RN INT06110
- RL3 = RL*RL*RL INT06120
- C = 4.*RL3*FF/RN INT06130
- CC = 0.0 INT06140
- A = 1.0 INT06150
- IF(IABS(MODE)-2)202,200,201 INT06160
- 200 CC = 0.5*PI*DR*DR*DR*F(1) INT06170
- A = .125*RN*RN*RN/(RL3*RL3) INT06180
- GOTO 202 INT06190
- 201 CC = F0 INT06200
- 202 IF(MODE.GT.0)GOTO 203 INT06210
- C = 1.0/C INT06220
- CC = A*CC INT06230
- 203 CONTINUE INT06240
- C COPY INPUT VECTOR INTO WORKING STORAGE INT06250
- DO 204 I=1,N INT06260
- 204 W(ND+I+1) = C*FF*F(I)*I INT06270
- C Eric Starts debugging here
- C RETURN
- C END
- CALL RSA(M) INT06280
- C COPY INTO OUTPUT VECTOR
- C N eq 2048
- C DO 205 I=1,N
- DO 205 I=1,1
- 205 T(I) = W(NR+I+1)/FLOAT(I)+CC INT06310
- RETURN INT06320
- END INT06330
- C**************************** SETFT **********************************
- SUBROUTINE SETFT(MM,M) INT06340
- C PREPARE FIXED DATA FOR FAST FOURIER TRANSFORM INT06350
- COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES INT06360
- COMMON /FWORK/ W(12288) INT06370
- IF(MM.LT.0)RETURN INT06380
- IF(MM.EQ.M)RETURN INT06390
- M = MM INT06400
- N = 2**M INT06410
- N2 = N/2 INT06420
- F = 1.0/FLOAT(N) INT06430
- DA = 6.2831853071796D0*F INT06440
- RTTWO = SQRT(2.0) INT06450
- F = SQRT(F) INT06460
- N1 = (N2-1)/2 INT06470
- NC = 0 INT06480
- NS = N1 INT06490
- NDATA = N1+N1 INT06500
- IF(N1.LE.0)RETURN INT06510
- DO 10 I=1,N1 INT06520
- A = FLOAT(I)*DA INT06530
- W(I) = COS(A) INT06540
- 10 W(N1+I) = SIN(A) INT06550
- RETURN INT06560
- END
- C************************** RSA *************************************
- SUBROUTINE RSA(M) INT06580
- C REAL ODD ANALYSIS OR SYNTHESIS INT06590
- COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES INT06600
- COMMON /FWORK/ W(12288) INT06610
- NQ = N2 INT06620
- NIB = NDATA +1 INT06630
- NOB = NDATA+N2+2 INT06640
- DO 200 L = 1,M INT06650
- NIM = NIB+NQ INT06660
- NIE = NIM+NQ INT06670
- NQH = NQ/2 INT06680
- NQL = (NQ-1)/2 INT06690
- NOM = NOB+NQL+1 INT06700
- C TPRIME = 0 INT06710
- IF(NQL.EQ.0)GOTO 20 INT06720
- DO 10 IT = 1,NQL INT06730
- W(NOB+IT)=W(NIB+IT)-W(NIM-IT) INT06740
- 10 W(NOM+IT) = W(NIB+IT)+W(NIM-IT) INT06750
- 20 IF(NQH.NE.NQL) W(NOM) = 2.0*W(NIB+NQH)
- C ?? 0.LT.TPRIME.LT.NP/2. INT06770
- NO1 = NQ INT06780
- NO2 = N2-NQ INT06790
- IF(NO1-NO2)40,80,120 INT06800
- 40 NI=NO1+NO1 INT06810
- W(NOB+NO1) = W(NIB+NI)+W(NIM+NI) INT06820
- W(NOB+NO2) = -W(NIB+NI)+W(NIM+NI) INT06830
- IF(NQL.EQ.0)GOTO 60 INT06840
- CC = W(NC+NO1) INT06850
- SS = W(NS+NO1) INT06860
- DO 50 IT = 1,NQL INT06870
- RE = CC*W(NIE+NI-IT)-SS*W(NIM+NI-IT) INT06880
- AI = SS*W(NIE+NI-IT)+CC*W(NIM+NI-IT) INT06890
- W(NOM+NO1+IT) = W(NIM+NI+IT)-RE INT06900
- W(NOM+NO2+IT) = W(NIM+NI+IT)+RE INT06910
- W(NOB+NO1+IT) = W(NIB+NI+IT)+AI INT06920
- 50 W(NOB+NO2+IT) = -W(NIB+NI+IT)+AI INT06930
- 60 IF(NQH.EQ.NQL)GOTO 70 INT06940
- NR = NO1/2 INT06950
- W(NOM+NO1) = 2.0*(W(NS+NR)*W(NIM+NI+NQH)+W(NC+NR)*W(NIB+NI+NQH)) INT06960
- W(NOM+NO2) = 2.0*(W(NC+NR)*W(NIM+NI+NQH)-W(NS+NR)*W(NIB+NI+NQH)) INT06970
- 70 NO1 = NO1+NQ INT06980
- NO2 = NO2-NQ INT06990
- IF(NO1-NO2)40,80,120 INT07000
- 80 W(NOB+NO1)=W(NIM) INT07010
- IF(NQL.EQ.0)GOTO 100 INT07020
- DO 90 IT = 1,NQL INT07030
- W(NOM+NO1+IT) = W(NIM+IT) INT07040
- 90 W(NOB+NO1+IT) = W(NIE-IT) INT07050
- 100 IF(NQH.NE.NQL)W(NOM+NO1) = RTTWO*W(NIM+NQH) INT07060
- 120 NT = NIB INT07070
- NIB = NOB INT07080
- NOB = NT INT07090
- NQ = NQH INT07100
- 200 CONTINUE INT07110
- NRES = NIB-1 INT07120
- RETURN INT07130
- END INT07140
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement