Appendix - Program Listings
The programs listed here were written using FORTRAN 77 and compiled using both Salford F77 and Acorn Desktop FORTRAN 77. The programs were written by both Dr. Trevor Crowley and myself and further developed over time.
Results for FindEE9 can be found in the results section and also in Capillary Condensation in a Porous Medium by Annmarie Moran, BSc final year project, University of Salford, 1999.
Full sources are available from here
Liquid Bridge - FindEE9
C---- MY VERSION OF PAUL'S FINDEE8.
C--- USING MY STANDARD OUTPUT ROUTINES
C---- MODIFIED FINDEE7 BY REPLACING BRIDGE()
C---- BY LATEST VERSION, ELIMINATING AREA/VOLUME BUG
C ---- VERSION 9.1. THIS VERSION CALLS SMAX BEFORE ENTERING S.
C OPTIONS(CHECK)
PROGRAM MAIN
C---- VERSION 9.1
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/BRIDGEPAR/FAL,FAW,FVT,FVCAP,FVL,SVT
COMMON/ENF/E,F
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,SNU,
1 SN2U,DNU,DN2U
COMMON/CURVE/H,KR0,KZ0,KRC,KZC
COMMON/RZ/RR,ZZ
DOUBLE PRECISION RR,ZZ
DOUBLE PRECISION H,KR0,KZ0,KRC,KZC
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION E,F
DOUBLE PRECISION FAL,FAW,FVT,FVCAP,FVL,SVT
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION SM ,RCMAX,DRC
DOUBLE PRECISION FVLVEC,FAWVEC,FALVEC,RCVEC
INTEGER LL,NRC,IS,II
CHARACTER*2 SYSTEM,YN
DOUBLE PRECISION PI
DIMENSION FVLVEC(1000),FAWVEC(1000),FALVEC(1000),SM(1000)
DIMENSION rcvec(1000)
CHARACTER*8 FNAME
CHARACTER*12 NAME
CHARACTER*4 EXT
C
C RTS,S & RC are both user inputted.
C
C Test value : RC=2, RTS=1 and S=0 should give EEO as 1.57398
C (via Mathematica 7)
C AL = 1.25603
C VT = 0.342391
C(via drop1.for)
PI = 3.14159D0
S=0.D0
WRITE (*,'(''*******************************************'')')
WRITE (*,*)
WRITE (*,'('' FRACTIONAL AREAS FOR LIQUID BRIDGES'')')
WRITE (*,*)
WRITE (*,'('' VERSION 9.1'')')
WRITE (*,*)
WRITE (*,'('' BY PAUL F. JOHNSON, 6/98 TO 3/99'')')
WRITE (*,*)
WRITE (*,'(''*******************************************'')')
WRITE (*,*)
WRITE (*,'('' PLEASE ENTER RTS : '',$)')
READ *,RTS
10 WRITE (*,'('' PLEASE ENTER RCMAX (-VE TO STOP),NRC: '',$)')
READ *,RCMAX,NRC
IF(RCMAX.LE.0.)STOP
DRC = RCMAX/NRC
II=1
IS=1
C PRINT *,'SMAX = ',SMAX
C PRINT *,'THETA1 = ',THETA1
C WRITE (*,*)
C WRITE (*,'('' FOR YOUR VALUE OF Rts, THE MAXIMUM SEPARATION'')')
C WRITE (*,'('' FOR THE SPHERE CAN BE '',$)')
C WRITE (*,*)SMAX
20 WRITE (*,'('' PLEASE ENTER S : '',$)')
READ *,S
IF (S.GT.SMAX) THEN
WRITE (*,*)
WRITE (*,'('' YOUR SMAX FOR THIS VALUE IS '',SMAX)')
GOTO 20
ENDIF
WRITE (*,'('' PLEASE ENTER SIN(BETA) : '',$)')
READ*,SB
WRITE (*,'('' SINGLE SPHERE OR 12 NEIGHBOURS (S OR N) '',$)')
READ (*,'(A)') SYSTEM
is=1
IF (SYSTEM.EQ.'N'.OR.SYSTEM.EQ.'n') IS=12
WRITE (*,'('' DO YOU WANT FVL TO BE MULTIPLIED BY 10 (Y/N)'',$)')
READ (*,'(A)') YN
ii = 1
IF (YN.EQ.'Y'.OR.YN.EQ.'y') II=10
DO 100 LL=1, NRC
RCAP = LL*DRC
CALL FINDSMAX0()
IF (S.GT.SMAX) THEN
WRITE (*,'('' S IS TOO LARGE FOR A BRIDGE TO FORM'')')
WRITE (*,*)
GOTO 20
ENDIF
CALL BRIDGE()
rcvec(ll) = rcap
SM(LL)=SMAX
FAW= FAW*IS
FAL= FAL*IS
FVL= FVL*II*IS
FAWVEC(LL)= FAW
FALVEC(LL)= FAL
FVLVEC(LL)= FVL
100 CONTINUE
WRITE (*,'('' PLEASE ENTER THE FILENAME (8 CHARS) '',$)')
READ (*,'(A)') FNAME
EXT='.CSV'
name=FNAME//EXT
CALL SAVER(RCVEC,FAWVEC,FALVEC,FVLVEC,SM,NRC,NAME)
PRINT *,'Rcap = ',rcap
PRINT *,'RN = ',RN0
PRINT *,'ZCS = ',ZCS
PRINT *,'RCS = ',RCS
PRINT *,'H = ',-1./RCAP
PRINT *,'FAL, AL = ', FAL, 4*PI*RTS**2*FAL
PRINT *,'FAW, AW = ', FAW, 4*PI*RTS**2*FAW
PRINT *,'FVT, VT = ', FVT, (4./3.)*PI*RTS**3*FVT
PRINT *,'FVCAP, VCAP = ', FVCAP, (4./3.)*PI*RTS**3*FVCAP
PRINT *,'FVL, VL = ', FVL, (4./3.)*PI*RTS**3*FVL
GOTO 10
END
C**************************************************************
c----- ROUTINES TO OBTAIN NECK PROFILE OF BRIDGE
C AND CURVATURE
C**************************************************************
DOUBLE PRECISION FUNCTION RN(ZZ)
C------ TO CALCULATE THE NECK PROFILE RN AS A FUNCTION OF Z
C----- HAVING OBTAINED A AND ES BY CALL TO BRIDGE
C----- TYPICALLY FOR 0 < Z < ZCS
C------ (AS IN GENERATING PROFILE IMPLICITLY GENERATE LOCAL
C----- VALUES OF ELLIPTIC FNS + INTERGRALS E.G. SNU = SIN(PHI(Z))
C---- E,F ETC IN COMMON /PARM/ AND /ENF/
C------ IF REQ'D COULD EASILY OBTAIN AL(Z),VT(Z), S(Z),PHI(Z) ETC.
C--- AFTER CALL TO RN(Z) BY EVALUATING THEM AS IN BRIDGE() )
COMMON/CNUS/CNUS,KS
COMMON/RZ/R,Z0
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION ZZ,Z,R,Z0,ZFN,PHIZ
DOUBLE PRECISION CUP,CLOW, XACC,RTBIS
EXTERNAL ZFN
IF(ZZ.GT.0.D0)THEN
Z = ZZ
ELSE
Z = -ZZ
ENDIF
XACC=1.D-7
CLOW= CNUS
CUP= 1.D0 - 1.D-7
PHIZ=RTBIS(ZFN,Z,CLOW,CUP,XACC)
RN = R
RETURN
END
DOUBLE PRECISION FUNCTION ZFN(CNUU)
C------ TO CALCULATE ZFN(SNU)
C------ TO GENERATE PROFILE
C----- BY FINDING ROOT TO Z = ZFN(SNU)
C----- HAVING A AND ES OBTAINED BY CALL TO BRIDGE
COMMON/RZ/R,Z
DOUBLE PRECISION CNUU,R,Z
CALL RANDZ(CNUU)
ZFN = Z
RETURN
END
SUBROUTINE RANDZ(CNUU)
COMMON/RZ/R,Z
C---- ROUTINE TO CALCULATE PARAMETRIC PROFILE
C----- R,Z SIMULTANEOUSLY AS A FUNCTION OF SNU = CNUU
C----- GIVEN A AND ES OBTAINED BY CALL TO BRIDGE
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION R,Z,ADK
DOUBLE PRECISION CNUS,KS,CNUU
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION L
DOUBLE PRECISION E,F
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION TNU
CALL ENF1(CNUU,KS)
ADK = A/K
c-- r(phi) = a ( k-1dn u - cn u )
R = ADK *(DNU - K*CNU)
C--- L = ARC LENGTH OF HYPERBOLA
c l(phi) = a k-1 (dn u tn u - E(phi,k) + k'2 F(phi,k))
TNU = SNU/CNU
L = ADK*(KDASH2*F-E+DNU*TNU)
c z(phi) = l(phi) - r(phi) tn u
Z = L - R*TNU
RETURN
END
DOUBLE PRECISION FUNCTION RSPH(ZZ)
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
C------ TO CALCULATE THE SPHERE PROFILE RS AS A FUNCTION OF Z
C----- TYPICALLY FOR S < Z < 2*RTS + S
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION ZZ,Z,ZSPH
IF(ZZ.GT.0.D0)THEN
Z = ZZ
ELSE
Z = -ZZ
ENDIF
IF(Z.LE.S.OR.Z.GE.RTS+RTS+S)THEN
RSPH = 0.D0
RETURN
ELSE
ZSPH = RTS + S - Z
RSPH = DSQRT(RTS*RTS- ZSPH*ZSPH)
ENDIF
RETURN
END
SUBROUTINE CURVEZ(KR,KZ,Z)
C------ TO CALCULATE THE RADIAL AND Z CURVATURES
C------- A FUNCTION OF Z - MAINLY TO BE USED FOR CHECKING
C----- HAVING OBTAINED A AND ES BY CALL TO BRIDGE
C----- TYPICALLY FOR 0 < Z < ZCS
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION KR,KZ
DOUBLE PRECISION Z,R,H, RN
C DOUBLE PRECISION DZ,SNU0, DSNU
c DOUBLE PRECISION KZ1
C---- H = -1/rcap = (1/2)(KZ+ KR ),
H = -1.D0/RCAP
R = RN(Z)
C----
KR = CNU/R
C----
KZ = 2.D0*H - KR
C---- KZ1 = -d SNU/dz for checking
c SNU0 = SNU
c DZ = 1.D-4*ZCS
c R = RN(Z+ DZ)
c DSNU = SNU - SNU0
c KZ1 = - DSNU/DZ
RETURN
END
C**************************************************************
c----- ROUTINES TO OBTAIN BASIC BRIDGE PARAMETERS
C----- AND GEOMETRY AT CONTACT OF BRIDGE WITH SPHERE
C**************************************************************
SUBROUTINE BRIDGE()
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/BRIDGEPAR/FAL,FAW,FVT,FVCAP,FVL,SVT
COMMON/CURVE/H,KR0,KZ0,KRC,KZC
DOUBLE PRECISION H,KR0,KZ0,KRC,KZC
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,
1 SNU,SN2U,DNU,DN2U
C------ TO CALCULATE SHAPE PARAMETERS OF BRIDGE
C------- IF S > 0, NEED TO HAVE CALLED FINDSMAX FIRST
C------ GIVEN RTS, S, RCAP IN COMMON /ARTS/
c --- where rts = radius of tangent sphere (beware = RCAP in scattering papers)
c----- s = (1/2) separation of surfaces of spheres (del0 = 2s )
c------ rcap = capillary radius = 1/H, H = mean curvature
C------ RETURNS (A)IN COMMON/ARTS/
C------ A,RN0,RCS,ZCS
c------ where: a = semi-axis of generating hyperbola = (1/2)rcap
c------ rn0 = neck radius
c------ rcs = contact radius
c------ zcs = z for contact (z= 0 at neck)
C------(B)IN COMMON/BRIDGEPAR/
C------ (i)FRACTIONAL AREAS FAL, FAW
C------ OF LIQUID INTERFACE, AL, AND WET SOLID, AW,
C----- SCALED TO AREA OF SPHERE,
C------ (ii)FRACTIONAL VOLUMES FVT, FVCAP, FVL
C------ OF TOTAL BRIDGE, VT, AND SOLID ENDCAP, VCAP,
C----- AND LIQUID VOLUME, VL,
C------ SCALED TO VOLUME OF SPHERE
C------ (iii)VIRIAL FUNCTION OF TOTAL FLUID BRIDGE, SVT,
C------ SCALED TO (VOLUME OF SPHERE/RTS)
C------ (i.e.excluding endcaps,for unit surface tension)
c-------- St = (1/3)(3a-1Vt + 2Al ), svt = st / (Vsphere/rts)
c-------- note: dp = 2 g H = g a-1
C--------- (C)IN COMMON/CURVE/H,KR0,KZ0,KRC,KZC
C----- MEAN CURVATURE AND CURVATURES AT
C----- NECK (Z=0) AND CONTACT (Z=RCS)
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION E,F
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION FAL,FAW,SVT,FVT,FVCAP,FVL
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION SNUM1,SINPHIBC,COSPHIBC,ARTS,SRTS,Spbc1
IF(S.GT.1.D-6)THEN
IF(SMAX-S.LT.1.D-6)THEN
C--- NO BRIDGE FORMED
KS = 1.D0
RN0= 0.D0
RCS= 0.D0
ZCS= 0.D0
FAL = 0.D0
FAW = 0.D0
SVT = 0.D0
FVT = 0.D0
FVCAP = 0.D0
FVL = 0.D0
RETURN
ENDIF
ENDIF
C--- BRIDGE FORMED
CALL FINDKS()
ARTS = A/RTS
SRTS = S/RTS
SNUM1 = 1.D0-SNU
IF(SB.LE.0.D0)THEN
ZCS = S + RTS *SNUM1
RCS = RTS*CNU
SINPHIBC = SNU
COSPHIBC = CNU
ELSE
SINPHIBC = SNU*CB+CNU*SB
COSPHIBC = CNU*CB-SNU*SB
ZCS = S + RTS *(1.D0 - SINPHIBC)
RCS = RTS*COSPHIBC
ENDIF
Spbc1 = 1.D0-SINPHIBC
RN0 = A*((1.D0/KS)-1.D0)
c-- Al= 2pi a2 k-1 (2 E(phis,k) - k'2 F(phiss,k) - 2 k sn u )
FAL = 2.D0*E - KDASH2* F - 2.D0*K *SNU
FAL= (0.5D0* (ARTS)**2 /K) *FAL
c FAW = 0.5D0*SNUM1
FAW = 0.5D0*Spbc1
c--- evaluate Vt using virial theorem where
C St = (1/3)(3a-1Vt + 2Al )
C = (1/3)(pi a (k'/k)2 z(phis) + 2pi snu r2(phis))
C SVT = RTS*ST/VSPHERE
C SVT = 0.25D0*(ARTS*(KDASH2/K2)*(SRTS+SNUM1)+ 2.0*SNU*CN2U)
SVT = 0.25D0*(ARTS*(KDASH2/K2)*(SRTS+SPBC1)
1 + 2.0*SNU*COSPHIBC*COSPHIBC)
FVT = ARTS *(SVT -(FAL+FAL))
C FVCAP = 0.25D0* SNUM1*SNUM1*(2.D0 + SNU)
FVCAP = 0.25D0* Spbc1*Spbc1*(2.D0 + SINPHIBC)
FVL = FVT-FVCAP
C---- CALCULATE CURVATURES AT CONTACT AND ZERO
C---- H = -1/rcap = (1/2)(KZ+ KR ),
H = -1.D0/RCAP
KR0 = 1.D0/RN0
KZ0 = 2.D0*H - KR0
IF(SB.LE.0.D0)THEN
KRC = 1.D0/RTS
ELSE
KRC = 1.D0/(RTS*(CB-SB*SNU/CNU))
ENDIF
KZC = 2.D0*H - KRC
RETURN
END
SUBROUTINE FINDKS()
C----- GIVEN RTS,S,SB,CB,RCAP IN COMMON/ARTS
C------- IF S > 0, NEED TO HAVE CALLED FINDSMAX FIRST
C---- RETURNS CNUS,KS IN COMMON/CNUS/
C------ WHERE KS IS ECCENTRICITY OF REQUIRED NODOID
C------AND CNUS IS VALUE OF SNU AT CONTACT
C---- AND CREATES ELLIPTIC PARAMS IN COMMON/PARM/
C----- AND E,F IN COMMON/ENF/
C----- FOR USE IN SUBROUTINE BRIDGE ETC.
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/CNUS/CNUS,KS
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION KS,CNUS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION SFN1,XACC,DEPS,RTBIS
EXTERNAL SFN1
XACC=1.D-10
DEPS =1.D-6
CNUS=RTBIS(SFN1,S,CNUMAX,1.D0-DEPS,XACC)
KS = K
RETURN
END
C**************************************************************
c----- ROUTINES FOR SMAX
C**************************************************************
SUBROUTINE FINDSMAX0()
C----- GIVEN RTS,RCAP IN COMMON/ARTS
C---- RETURNS CNUMAX, SMAX IN COMMON/SMAX/
C----- AND CALCULATES A
C---- THIS VERSION USES MAXTRI ON SFN1
C--- FOR ROBUSTNESS/CHECKING
C----- SMAX USED BY BRIDGE TO CHECK IF BRIDGE FORMED
c---- CNUMAX USED BY FINDKS -> LIMITS FOR ROOT SEARCH
C------ USES SMAXFN
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION SFN1,MAXTRI
DOUBLE PRECISION XACC,DEPS
EXTERNAL SFN1
A = 0.5D0 * RCAP
DEPS =1.D-6
XACC=1.D-10
C------ CAN IMPROVE ON LIMITS FOR ROOT LATER
C---- FOR FINITE CONTACT ANGLE
CNUMAX =MAXTRI(SFN1,DEPS,1.D0-DEPS,XACC)
SMAX= SFN1(CNUMAX)
RETURN
END
DOUBLE PRECISION FUNCTION SFN1(CNUU)
C------ TO CALCULATE SFN(CNU) FOR ROOT FINDING
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION KFN,CNUU,L,KK
DOUBLE PRECISION E,F
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
KK = KFN(CNUU)
IF(KK.GE.1.D0)THEN
C---- IF KFN OR INVALID SET SFN1 VERY NEGATIVE, FOR MAX FINDING
SFN1 = - RTS
ELSE
CALL ENF1(CNUU,KK)
L=(A/K)*(KDASH2*F-E+DNU*SNU/CNU)
SFN1=L-RTS*(1.D0 - SB/CNUU)
ENDIF
RETURN
END
C**************************************************************
c----- ROUTINE FOR K AS FN OF CNU
c---- analogous to ee fn of theta
C**************************************************************
DOUBLE PRECISION FUNCTION KFN(CNU)
C----- BASIC CONTACT FUNCTION TO GIVE K AS FN OF CNU
C---- USED IN SFN1 FOR FINDKS
C---- FOR CONTACT ANGLE BETA, SB = SIGN BETA
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RS,ZS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RS,ZS
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,FF,RTSDA
C----- NOTE REDUNDANCY HERE AS TYPICALLY
C----- WHEN CALL KFN WILL ALSO CALL ENF1
C----- SUCH THAT CN2U,SN2U,SNU CALCULATED TWICE
C----- TO AVOID THIS WOULD BE COMPLICATED
C----- AND WOULD MAKE ENF1 LESS SIMPLE FOR USE BY OTHER ROUTINES
C----- COULD MODIFY LATER BY SPLITTING SNU ETC, FROM K AND E+ F EVALN
C----- WHICH COULD BE USED IN SFN1
CB =DSQRT(1.D0 -SB*SB)
CN2U = CNU*CNU
SN2U = 1.D0-CN2U
SNU = DSQRT(SN2U)
RTSDA= RTS/A
IF(SB.LE.0.D0)THEN
KFN = DSQRT(1.D0/((1.D0+RTSDA)**2*CN2U+SN2U))
ELSE
FF = (RTSDA*CB+1.D0)*CNU - RTSDA*SB*SNU
KFN = DSQRT(1.D0/(FF**2+SNU*SNU))
ENDIF
RETURN
END
C**************************************************************
c----- ROUTINES FOR CNUMAXMAX
C**************************************************************
SUBROUTINE FINDCNUMAXMAX()
C----- CALCULATES CNUMAXMAX, SMAXMAX
C------- CORRESPONDING TO CNUMAX, SMAX FOR RCAP = INFINITE
C---- ALSO CALCULATES CB = COS(BETA)
C---- CNUMAXMAX USED BY FINDSMAX, FINDKS
C---- -> LIMITS FOR ROOT SEARCH
C------ USES CNUMAXMAXFN
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/CNUMAXMAX/CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CNU,SNU,COSPHIB,SIN2PHIB,SINPHIB
DOUBLE PRECISION CNUMAXMAXFN
DOUBLE PRECISION XACC,RTBIS
EXTERNAL CNUMAXMAXFN
IF(SB.LE.0.D0)THEN
CB = 1.D0
CNUMAXMAX = 0.552434135D0
SMAXMAX = 0.199678650D0
ELSE
CB = DSQRT(1.D0 - SB*SB)
XACC=1.D-10
C------ CAN IMPROVE ON LIMITS FOR ROOT LATER
C------ USE BISECTION AS SECANT DOES NOT WORK STABLELY
CNUMAXMAX=RTBIS(CNUMAXMAXFN,0.D0,0.552434135D0,1.D0-1.D-7,XACC)
C--- NOW CALCULATE SMAXMAX
C--- COSPHIB = COS(PHI+BETA) - CALCULATED IN CNUMAXMAXFN
C------ CNU, SNU - CALCULATED IN CNUMAXMAXFN
C--- SIN2PHIB = COS(2*PHI+BETA) - CALCULATED IN CNUMAXMAXFN
C--- SINPHIB = SIN(PHI+BETA)
SINPHIB = SNU*CB + CNU*SB
SMAXMAX = 2.D0*CNU* COSPHIB* COSPHIB/SIN2PHIB
SMAXMAX = RTS*(SMAXMAX + SINPHIB -1.D0)
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION CNUMAXMAXFN(CNUU)
C------ TO CALCULATE CNUMAXMAXFN(CNU) FOR FINDING CNUMAXMAX
C---- FOR SNU = CNUU, K = KFN(CNUU)
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/CNUMAXMAX/CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CN2U,SN2U,CNUU
DOUBLE PRECISION DACOSH
CNU = CNUU
CN2U = CNU*CNU
SN2U = 1.D0 -CN2U
SNU = DSQRT(SN2U)
C--- COSPHIB = COS(PHI+BETA)
COSPHIB = CNU*CB - SNU*SB
C--- SIN2PHIB = COS(2*PHI+BETA)
SIN2PHIB = 2.D0*SNU*CNU*CB + (CN2U-SN2U)*SB
C CNUMAXMAXFN = DCOSH(2.D0*COSPHIB/SIN2PHIB)-1.D0/CNU
CNUMAXMAXFN = (2.D0*COSPHIB/SIN2PHIB)-DACOSH(1.D0/CNU)
C PRINT *,CNUMAXMAXFN
RETURN
END
C**************************************************************
c----- ROUTINES FOR ELLIPTIC PARAMETERS AND INTEGRALS
C**************************************************************
SUBROUTINE ENF1(CNUU,KK)
c-------- main private routine - together with parmgen
c---- given K=KK and SNU = CNUU
c----- for some applications generates more info than required
c----- but efficient - so don't modify
COMMON/ENF/E,F
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
C----- RETURNS E AND F IN COMMON
C---- GENERATES ELLIPTIC FUNCTIONS etc. BY CALL TO
C------ PARMGEN
C----- and GENERATES E AND F USING CARLSONS RD, RF FUNCTIONS
DOUBLE PRECISION ERRTOL, RF, RD, F, E
DOUBLE PRECISION CNUU,KK
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
INTEGER IERR
CALL PARMGEN(CNUU,KK)
ERRTOL=0.001
F = SNU*RF(CN2U,DN2U,1.D0,ERRTOL,IERR)
E =F-(1.D0/3.D0)*K2*SNU*SN2U*RD(CN2U,DN2U,1.D0,ERRTOL,IERR)
RETURN
END
SUBROUTINE PARMGEN(CNUU,KK)
c--- subroutine to generate elliptic modulus, k, complementary modulus, k'
c--- and Jacobian elliptic functions together with their squares,
c---- given K=KK and CNU = CNUU
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,SNU,
1 SN2U,DNU,DN2U
DOUBLE PRECISION CNUU,KK
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
CNU = CNUU
K = KK
K2 = K * K
KDASH2 = 1.D0- K2
KDASH = DSQRT(KDASH2)
CN2U = CNU*CNU
SN2U = 1.D0 -CN2U
SNU = DSQRT(SN2U)
DN2U = 1.D0 -K2*SN2U
DNU = DSQRT(DN2U)
RETURN
END
C**************************************************************
c----- STANDARD ROUTINES
C**************************************************************
DOUBLE PRECISION FUNCTION MAXTRI(FN,X1,X2,XACC)
C---- TO FIND MAX BETWEEEN X1 AND X2
C---- BY TRISECTION OR SOMETHING SIMILAR
C--- SHOULD BE ROBUSt
c---- note: accuracy in xx = maxtri limited by accuracy in FN about max
c----- would obtain better accuracy in xx by finding root of derivative fn.
c--- however it is usually accuracy in FMAX = FN(MAXTRI) which is important
c--- could modify routine (and RTSEC/RTBIS) to return on accuracy of fn
c--- rather than xx
DOUBLE PRECISION X1,X2,FN,XACC,FMAX
DOUBLE PRECISION XOLD(3),FOLD(3),XNEW(5),FNEW(5),DX,XMID
INTEGER MAXIT
INTEGER LL,L1,L2,LMAX
PARAMETER (MAXIT=100)
DX = 0.5D0*(X2-X1)
XOLD(1) = X1
XMID = X1 + DX
XOLD(2) = XMID
XOLD(3) = X2
FOLD(1) = FN(X1)
FOLD(2) = FN(XMID)
FOLD(3) = FN(X2)
DO 100 LL =1,MAXIT
C--- CONSTRUCT 5 VALS BY BISECTION OF TWO INTERVALS
C--- FIRST COPY OLD TO POSN'S 1,3,5
DO 101 L1 = 1,3
L2 = L1+L1-1
XNEW(L2) = XOLD(L1)
FNEW(L2) = FOLD(L1)
101 CONTINUE
C--- NOW FILL POSN'S 2 AND 4
DX= 0.5D0*DX
DO 102 L1 = 1,2
L2 = L1 +L1
XMID=XOLD(L1) + DX
XNEW(L2) = XMID
FNEW(L2) = FN(XMID)
102 CONTINUE
C--- NOW FIND POSN. OF MAXIMIUM - MUST BE 2,3 OR 4
LMAX = 2
FMAX = FNEW(LMAX)
DO 103 L1 = 3,4
IF(FNEW(L1).GT.FMAX)THEN
LMAX = L1
FMAX = FNEW(LMAX)
ENDIF
103 CONTINUE
C----- FILL OLD WITH POINTS LMAX +/-1
DO 104 L1 = 1,3
L2 = LMAX + L1 -2
XOLD(L1) = XNEW(L2)
FOLD(L1) = FNEW(L2)
104 CONTINUE
C PRINT*, FMAX, XOLD(2),DX
IF (DX.LT.XACC)GOTO 10
100 CONTINUE
C PAUSE 'MAXTRI EXCEED MAX. ITERATIONS'
10 MAXTRI = XOLD(2)
RETURN
END
DOUBLE PRECISION FUNCTION RTBIS(FUNC, VAL,X1,X2,XACC)
C---- MODIFIED TO FIND ROOT FUNV = VAL
DOUBLE PRECISION X1,X2,FUNC,VAL,XACC
DOUBLE PRECISION FMID,F,XMID,DX
INTEGER JMAX,J
PARAMETER (JMAX=100)
FMID=FUNC(X2)-VAL
F=FUNC(X1)-VAL
IF (F.LT.0.D0) THEN
RTBIS=X1
DX = X2-X1
ELSE
RTBIS=X2
DX = X1-X2
ENDIF
DO 11 J=1,JMAX
DX= 0.5D0*DX
XMID=RTBIS + DX
FMID=FUNC(XMID)-VAL
C PRINT *,'HELLO',FMID
C READ*, LL
IF (FMID.LE.0.D0)RTBIS = XMID
IF (DABS(DX).LT.XACC) RETURN
11 CONTINUE
C PAUSE 'RTBIS EXCEED MAX. ITERATIONS'
END
DOUBLE PRECISION FUNCTION RTSEC(SFN1,S,TLOW,TUP,XACC)
C---- MODIFIED TO FIND ROOT SFN1 = S
C---- PREVIOUSLY ASSUMED SVAL = 0
DOUBLE PRECISION TLOW,TUP,XACC,XL,SFN1,S
DOUBLE PRECISION FL,F,SWAP,DX
INTEGER MAXIT,H
PARAMETER (MAXIT=30)
FL=SFN1(TLOW)-S
F=SFN1(TUP)-S
IF (ABS(FL).LT.ABS(F)) THEN
RTSEC=TLOW
XL=TUP
SWAP=FL
FL=F
F=SWAP
ELSE
XL=TLOW
RTSEC=TUP
ENDIF
DO 11 H=1,MAXIT
DX=(XL-RTSEC)*F/(F-FL)
XL=RTSEC
FL=F
RTSEC=RTSEC+DX
F=SFN1(RTSEC) -S
IF (DABS(DX).LT.XACC) RETURN
11 CONTINUE
C PAUSE 'RTSEC EXCEED MAX. ITERATIONS'
END
DOUBLE PRECISION FUNCTION RD(X,Y,Z,ERRTOL,IERR)
INTEGER IERR
DOUBLE PRECISION C1,C2,C3,C4,EA,EB,EC,ED,EF,EPSLON,ERRTOL,LAMDA
DOUBLE PRECISION LOLIM,MU,POWER4,SIGMA,S1,S2,UPLIM,X,XN,XNDEV
DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
DATA LOLIM/1.0D-25/,UPLIM/4.5D+21/
IF (DMIN1(X,Y).LT.0.D0) GOTO 100
IF (DMIN1(X+Y,Z).LT.LOLIM) GOTO 100
IF (DMAX1(X,Y,Z).LE.UPLIM) GOTO 112
100 PRINT *,'INVALID ARGUMENTS PASSED TO RD'
IERR=1
GOTO 124
112 IERR=0
XN=X
YN=Y
ZN=Z
SIGMA=0.D0
POWER4=1.D0
116 MU=(XN+YN+3.D0*ZN)*0.2D0
XNDEV=(MU-XN)/MU
YNDEV=(MU-YN)/MU
ZNDEV=(MU-ZN)/MU
EPSLON=DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
IF (EPSLON.LT.ERRTOL) GOTO 120
XNROOT=DSQRT(XN)
YNROOT=DSQRT(YN)
ZNROOT=DSQRT(ZN)
LAMDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
SIGMA=SIGMA+POWER4/(ZNROOT*(ZN+LAMDA))
POWER4=POWER4*0.25D0
XN=(XN+LAMDA)*0.25D0
YN=(YN+LAMDA)*0.25D0
ZN=(ZN+LAMDA)*0.25D0
GOTO 116
120 C1=3.D0/14.D0
C2=1.D0/6.D0
C3=9.D0/22.D0
C4=3.D0/26.D0
EA=XNDEV*YNDEV
EB=ZNDEV*ZNDEV
EC=EA-EB
ED=EA-6.D0*EB
EF=ED+EC+EC
S1=ED*(-C1+0.25D0*C3*ED-1.5D0*C4*ZNDEV*EF)
S2=ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA))
RD=3.D0*SIGMA+POWER4*(1.D0+S1+S2)/(MU*DSQRT(MU))
124 RETURN
END
DOUBLE PRECISION FUNCTION RF(X,Y,Z,ERRTOL,IERR)
INTEGER IERR
DOUBLE PRECISION C1,C2,C3,E2,E3,EPSLON,ERRTOL,LAMBDA
DOUBLE PRECISION LOLIM,MU,S,UPLIM,X,XN,XNDEV,XNROOT
DOUBLE PRECISION Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
DATA LOLIM/1.5D-38/, UPLIM/ 3.0D+37/
IF (DMIN1(X,Y,Z).LT.0.D0) GOTO 100
IF (DMIN1(X+Y,X+Z,Y+Z).LT.LOLIM) GOTO 100
IF (DMAX1(X,Y,Z).LE.UPLIM) GOTO 112
100 PRINT *,'INVALID ARGUMENTS PASSED TO RF'
IERR=1
GOTO 124
112 IERR=0
XN=X
YN=Y
ZN=Z
116 MU=(XN+YN+ZN)/3.D0
XNDEV=2.D0-(MU+XN)/MU
YNDEV=2.D0-(MU+YN)/MU
ZNDEV=2.D0-(MU+ZN)/MU
EPSLON=DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
IF (EPSLON.LT.ERRTOL) GOTO 120
XNROOT=DSQRT(XN)
YNROOT=DSQRT(YN)
ZNROOT=DSQRT(ZN)
LAMBDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
XN=(XN+LAMBDA)*0.25D0
YN=(YN+LAMBDA)*0.25D0
ZN=(ZN+LAMBDA)*0.25D0
GOTO 116
120 C1=1.D0/24.D0
C2=3.D0/44.D0
C3=1.D0/14.D0
E2=XNDEV*YNDEV-ZNDEV*ZNDEV
E3=XNDEV*YNDEV*ZNDEV
S=1.D0+(C1*E2-0.1D0-C2*E3)*E2+C3*E3
RF=S/DSQRT(MU)
124 RETURN
END
DOUBLE PRECISION FUNCTION DACOSH(X)
DOUBLE PRECISION X
DACOSH = DLOG (X + DSQRT(X*X-1.D0))
RETURN
END
SUBROUTINE SAVER(RCVEC,FAWVEC,FALVEC,FVLVEC,SM,NRC,NAME)
INTEGER L,NRC
DOUBLE PRECISION FAWVEC(NRC),FALVEC(NRC),FVLVEC(NRC)
DOUBLE PRECISION SM(NRC),RCVEC(NRC)
CHARACTER*12 NAME
1 OPEN (4,FILE=NAME,ERR=999)
WRITE(4,11)
1 'RC','FAW','FAL','FVL','Smax'
11 FORMAT(5(1X,A,:,','))
DO 20 L=1,NRC
WRITE (4,12)
1 RCVEC(L),FAWVEC(L),FALVEC(L),FVLVEC(L),SM(L)
12 FORMAT(1P,5(G23.15,:,','))
20 CONTINUE
CLOSE (4)
RETURN
999 WRITE (*,'(''FILER ERROR.'')')
GOTO 1
END
Listing 1b - Bridge scattering - Britest3ba
PROGRAM MAIN
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/BRIDGEPAR/FAL,FAW,FVT,FVCAP,FVL,SVT
COMMON/ENF/E,F
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,SNU,
1 SN2U,DNU,DN2U
COMMON/CURVE/H,KR0,KZ0,KRC,KZC
COMMON/RZ/RR,ZZ
DOUBLE PRECISION RR,ZZ
DOUBLE PRECISION R,Z
DOUBLE PRECISION H,KR0,KZ0,KRC,KZC
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION E,F
DOUBLE PRECISION FAL,FAW,FVT,FVCAP,FVL,SVT
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION RN,RSPH,XACC,DEPS,RSP
DOUBLE PRECISION KR,KZ
DOUBLE PRECISION PI,NPT2,DZ,PHI
CHARACTER*10 NAME
INTEGER NPT,LL
REAL XVEC(1000),YVEC(1000)
DOUBLE PRECISION XX(100000),YY(100000),ZZV(100000)
INTRINSIC DSIN
PI = 4.D0*DATAN(1.D0)
DEPS = 1.D-3
XACC = 1.D-8
S=0.D0
WRITE (*,'('' PLEASE ENTER RADIUS OF SPHERE : '',$)')
READ *,RTS
5 WRITE (*,'('' PLEASE ENTER RCAP (-VE TO STOP): '',$)')
READ *,RCAP
A = 0.5D0*RCAP
IF(RCAP.LE.0.)STOP
WRITE (*,'('' PLEASE ENTER SIN(BETA) : '',$)')
READ *,SB
IF(SB.LT.0.)GOTO 5
CALL FINDCNUMAXMAX
CALL FINDSMAX0
WRITE (*,'('' PLEASE ENTER S : '',$)')
READ *,S
IF(S.LT.0.)GOTO 5
CALL BRIDGE()
CALL RANDZ(CNUS)
RSP = RSPH(Z)
NPT = 1000
NPT2=NPT/2
DZ = (S+RTS)/NPT2
DO 300 LL = 1,NPT
Z = (LL-NPT2-0.5)*DZ
XVEC(LL) = Z
300 CONTINUE
PRINT *,'STORING +VE BRIDGE PROFILE '
name='+VEBridge'
CALL FILEGEN(XVEC,YVEC,NPT,NAME)
DO 301 LL = 1,NPT
YVEC(LL) = -YVEC(LL)
301 CONTINUE
PRINT *,'STORING -VE BRIDGE PROFILE '
name='-VEBridge'
CALL FILEGEN(XVEC,YVEC,NPT,NAME)
DZ = (S+RTS)/NPT2
DO 320 LL = 1,NPT
Z = (LL-NPT2-0.5)*DZ
XVEC(LL) = Z
YVEC(LL) = RSPH(Z)
IF(Z.LE.0)YVEC(LL) = -YVEC(LL)
320 CONTINUE
PRINT *,'STORING +VE SPHERE PROFILE '
name='+VESphere'
CALL FILEGEN(XVEC,YVEC,NPT,NAME)
DO 321 LL = 1,NPT
YVEC(LL) = -YVEC(LL)
321 CONTINUE
PRINT *,'STORING -VE SPHERE PROFILE '
name='-VESphere'
CALL FILEGEN(XVEC,YVEC,NPT,NAME)
PRINT *,'CREATING 3D PROFILE - THIS MAY TAKE SOME TIME!'
NPT=100000
NPT2=NPT/2
PHI=0
DZ = (S+RTS)/NPT2
DO 400 L=1,NPT
ZZV(L) = (L-NPT2-0.5)*DZ
XX(L) = RN(ZZV(L))*DCOS(PHI)
IF(ABS(ZZV(L)).LE.ZCS) THEN
YY(L) = RN(ZZV(L))*DSIN(PHI)
ELSE
YY(L) = RSPH(ZZV(L))*DSIN(PHI)
ENDIF
PHI=PHI+(360./NPT)
IF (L.EQ.25000.) PRINT *,'25% DONE'
IF (L.EQ.50000.) PRINT *,'50% DONE'
IF (L.EQ.75000.) PRINT *,'75% DONE'
IF (PHI.GE.360.) GOTO 401
400 CONTINUE
401 PRINT *,'COMPLETE'
name='3DProfile'
CALL BIGSAVE(XX,YY,ZZV,NPT,NAME)
END
SUBROUTINE BIGSAVE(XX,YY,ZZV,NCH,NAME)
INTEGER NCH
CHARACTER*10 NAME
DOUBLE PRECISION XX(NCH),YY(NCH),ZZV(NCH)
INTEGER LL
OPEN(4,FILE=NAME,ERR=999)
WRITE(4,*)'X',',','Y',',','Z'
DO 1 LL=1,NCH
WRITE(4,*)XX(LL),',',YY(LL),',',ZZV(LL)
1 CONTINUE
CLOSE(4)
RETURN
999 WRITE(*,'('' ERROR IN TRYING THE FILE -TRY AGAIN !'')')
END
SUBROUTINE FILEGEN(X,Y,NCH,NAME)
INTEGER NCH
CHARACTER*10 NAME
REAL X(NCH),Y(NCH)
INTEGER LL
OPEN(4,FILE=NAME,ERR=999)
DO 1 LL=1,NCH
WRITE(4,*)X(LL),',',Y(LL)
1 CONTINUE
CLOSE(4)
RETURN
999 WRITE(*,'('' ERROR IN TRYING THE FILE -TRY AGAIN !'')')
END
DOUBLE PRECISION FUNCTION RN(ZZ)
COMMON/CNUS/CNUS,KS
COMMON/RZ/R,Z0
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION ZZ,Z,R,Z0,ZFN
DOUBLE PRECISION CUP,CLOW, XACC,RTBIS
DOUBLE PRECISION PHIZ
EXTERNAL ZFN
IF(ZZ.GT.0.D0)THEN
Z = ZZ
ELSE
Z = -ZZ
ENDIF
XACC=1.D-7
CLOW= CNUS
CUP= 1.D0 - 1.D-7
PHIZ=RTBIS(ZFN,Z,CLOW,CUP,XACC)
RN = R
RETURN
END
DOUBLE PRECISION FUNCTION ZFN(CNUU)
COMMON/RZ/R,Z
DOUBLE PRECISION CNUU,R,Z
CALL RANDZ(CNUU)
ZFN = Z
RETURN
END
SUBROUTINE RANDZ(CNUU)
COMMON/RZ/R,Z
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION R,Z,ADK
DOUBLE PRECISION CNUS,KS,CNUU
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION L
DOUBLE PRECISION E,F
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION TNU
CALL ENF1(CNUU,KS)
ADK = A/K
R = ADK *(DNU - K*CNU)
TNU = SNU/CNU
L = ADK*(KDASH2*F-E+DNU*TNU)
Z = L - R*TNU
RETURN
END
DOUBLE PRECISION FUNCTION RSPH(ZZ)
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION ZZ,Z,ZSPH
IF(ZZ.GT.0.D0)THEN
Z = ZZ
ELSE
Z = -ZZ
ENDIF
IF(Z.LE.S.OR.Z.GE.RTS+RTS+S)THEN
RSPH = 0.D0
RETURN
ELSE
ZSPH = RTS + S - Z
RSPH = DSQRT(RTS*RTS- ZSPH*ZSPH)
ENDIF
RETURN
END
SUBROUTINE CURVEZ(KR,KZ,Z)
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION KR,KZ
DOUBLE PRECISION Z,R,H, RN
H = -1.D0/RCAP
R = RN(Z)
KR = CNU/R
KZ = 2.D0*H - KR
RETURN
END
SUBROUTINE BRIDGE()
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/BRIDGEPAR/FAL,FAW,FVT,FVCAP,FVL,SVT
COMMON/CURVE/H,KR0,KZ0,KRC,KZC
DOUBLE PRECISION H,KR0,KZ0,KRC,KZC
COMMON/CNUS/CNUS,KS
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,
1 SNU,SN2U,DNU,DN2U
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION CNUS,KS
DOUBLE PRECISION E,F
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION FAL,FAW,SVT,FVT,FVCAP,FVL
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION SNUM1,SINPHIBC,COSPHIBC,ARTS,SRTS
IF(S.GT.1.D-6)THEN
IF(SMAX-S.LT.1.D-6)THEN
KS = 1.D0
RN0= 0.D0
RCS= 0.D0
ZCS= 0.D0
FAL = 0.D0
FAW = 0.D0
SVT = 0.D0
FVT = 0.D0
FVCAP = 0.D0
FVL = 0.D0
RETURN
ENDIF
ENDIF
CALL FINDKS()
ARTS = A/RTS
SRTS = S/RTS
SNUM1 = 1.D0-SNU
IF(SB.LE.0.D0)THEN
ZCS = S + RTS *SNUM1
RCS = RTS*CNU
ELSE
SINPHIBC = SNU*CB+CNU*SB
COSPHIBC = CNU*CB-SNU*SB
ZCS = S + RTS *(1.D0 - SINPHIBC)
RCS = RTS*COSPHIBC
ENDIF
RN0 = A*((1.D0/KS)-1.D0)
FAL = 2.D0*E - KDASH2* F - 2.D0*K *SNU
FAL= (0.5D0* (ARTS)**2 /K) *FAL
FAW = 0.5D0*SNUM1
SVT = 0.25D0*(ARTS*(KDASH2/K2)*(SRTS+SNUM1)+ 2.0*SNU*CN2U)
FVT = ARTS *(SVT -(FAL+FAL))
FVCAP = 0.25D0* SNUM1*SNUM1*(2.D0 + SNU)
FVL = FVT-FVCAP
H = -1.D0/RCAP
KR0 = 1.D0/RN0
KZ0 = 2.D0*H - KR0
IF(SB.LE.0.D0)THEN
KRC = 1.D0/RTS
ELSE
KRC = 1.D0/(RTS*(CB-SB*SNU/CNU))
ENDIF
KZC = 2.D0*H - KRC
RETURN
END
SUBROUTINE FINDKS()
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/CNUS/CNUS,KS
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION KS,CNU
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION SFN1,XACC,DEPS,RTBIS
EXTERNAL SFN1
XACC=1.D-10
DEPS =1.D-6
CNUS=RTBIS(SFN1,S,CNUMAX,1.D0-DEPS,XACC)
KS = K
RETURN
END
SUBROUTINE FINDSMAX0()
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
DOUBLE PRECISION SFN1,MAXTRI
DOUBLE PRECISION XACC,DEPS
EXTERNAL SFN1
A = 0.5D0 * RCAP
DEPS =1.D-6
XACC=1.D-10
CNUMAX =MAXTRI(SFN1,DEPS,1.D0-DEPS,XACC)
SMAX= SFN1(CNUMAX)
RETURN
END
DOUBLE PRECISION FUNCTION SFN1(CNUU)
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
COMMON/ENF/E,F
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION KFN,CNUU,L,KK
DOUBLE PRECISION E,F
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
KK = KFN(CNUU)
IF(KK.GE.1.D0)THEN
SFN1 = - RTS
ELSE
CALL ENF1(CNUU,KK)
L=(A/K)*(KDASH2*F-E+DNU*SNU/CNU)
SFN1=L-RTS*(1.D0 - SB/CNUU)
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION KFN(CNU)
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RS,ZS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RS,ZS
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,FF,RTSDA
CB =DSQRT(1.D0 -SB*SB)
CN2U = CNU*CNU
SN2U = 1.D0-CN2U
SNU = DSQRT(SN2U)
RTSDA= RTS/A
IF(SB.LE.0.D0)THEN
KFN = DSQRT(1.D0/((1.D0+RTSDA)**2*CN2U+SN2U))
ELSE
FF = (RTSDA*CB+1.D0)*CNU - RTSDA*SB*SNU
KFN = DSQRT(1.D0/(FF**2+SNU*SNU))
ENDIF
RETURN
END
SUBROUTINE FINDCNUMAXMAX()
COMMON/SMAX/CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION CNUMAX,CNUMAXMAX,SMAX,SMAXMAX
COMMON/CNUMAXMAX/CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CNU,SNU,COSPHIB,SIN2PHIB,SINPHIB
DOUBLE PRECISION CNUMAXMAXFN
DOUBLE PRECISION XACC,RTBIS
EXTERNAL CNUMAXMAXFN
IF(SB.LE.0.D0)THEN
CB = 1.D0
CNUMAXMAX = 0.552434135D0
SMAXMAX = 0.199678650D0
ELSE
CB = DSQRT(1.D0 - SB*SB)
XACC=1.D-10
CNUMAXMAX=RTBIS(CNUMAXMAXFN,0.D0,0.552434135D0,1.D0-1.D-7,XACC)
SINPHIB = SNU*CB + CNU*SB
SMAXMAX = 2.D0*CNU* COSPHIB* COSPHIB/SIN2PHIB
SMAXMAX = RTS*(SMAXMAX + SINPHIB -1.D0)
ENDIF
RETURN
END
DOUBLE PRECISION FUNCTION CNUMAXMAXFN(CNUU)
COMMON/ARTS/A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
DOUBLE PRECISION A,RTS,S,SB,CB,RCAP,RN0,RCS,ZCS
COMMON/CNUMAXMAX/CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CNU,SNU,COSPHIB,SIN2PHIB
DOUBLE PRECISION CN2U,SN2U,CNUU
DOUBLE PRECISION DACOSH
CNU = CNUU
CN2U = CNU*CNU
SN2U = 1.D0 -CN2U
SNU = DSQRT(SN2U)
COSPHIB = CNU*CB - SNU*SB
SIN2PHIB = 2.D0*SNU*CNU*CB + (CN2U-SN2U)*SB
CNUMAXMAXFN = (2.D0*COSPHIB/SIN2PHIB)-DACOSH(1.D0/CNU)
RETURN
END
SUBROUTINE ENF1(CNUU,KK)
COMMON/ENF/E,F
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,
1 CN2U,SNU,SN2U,DNU,DN2U
DOUBLE PRECISION ERRTOL, RF, RD, F, E
DOUBLE PRECISION CNUU,KK
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
INTEGER IERR
CALL PARMGEN(CNUU,KK)
ERRTOL=0.001
F = SNU*RF(CN2U,DN2U,1.D0,ERRTOL,IERR)
E =F-(1.D0/3.D0)*K2*SNU*SN2U*RD(CN2U,DN2U,1.D0,ERRTOL,IERR)
RETURN
END
SUBROUTINE PARMGEN(CNUU,KK)
COMMON/PARM/K,K2,KDASH,KDASH2,CNU,CN2U,SNU,
1 SN2U,DNU,DN2U
DOUBLE PRECISION CNUU,KK
DOUBLE PRECISION K,K2,KDASH,KDASH2
DOUBLE PRECISION CNU,CN2U,SNU,SN2U,DNU,DN2U
CNU = CNUU
K = KK
K2 = K * K
KDASH2 = 1.D0- K2
KDASH = DSQRT(KDASH2)
CN2U = CNU*CNU
SN2U = 1.D0 -CN2U
SNU = DSQRT(SN2U)
DN2U = 1.D0 -K2*SN2U
DNU = DSQRT(DN2U)
RETURN
END
DOUBLE PRECISION FUNCTION MAXTRI(FN,X1,X2,XACC)
DOUBLE PRECISION X1,X2,FN,XACC
DOUBLE PRECISION XOLD(3),FOLD(3),XNEW(5),FNEW(5),DX,XMID
INTEGER MAXIT,LL,L1,L2,LMAX
DOUBLE PRECISION FMAX
PARAMETER (MAXIT=100)
DX = 0.5D0*(X2-X1)
XOLD(1) = X1
XMID = X1 + DX
XOLD(2) = XMID
XOLD(3) = X2
FOLD(1) = FN(X1)
FOLD(2) = FN(XMID)
FOLD(3) = FN(X2)
DO 100 LL =1,MAXIT
DO 101 L1 = 1,3
L2 = L1+L1-1
XNEW(L2) = XOLD(L1)
FNEW(L2) = FOLD(L1)
101 CONTINUE
DX= 0.5D0*DX
DO 102 L1 = 1,2
L2 = L1 +L1
XMID=XOLD(L1) + DX
XNEW(L2) = XMID
FNEW(L2) = FN(XMID)
102 CONTINUE
LMAX = 2
FMAX = FNEW(LMAX)
DO 103 L1 = 3,4
IF(FNEW(L1).GT.FMAX)THEN
LMAX = L1
FMAX = FNEW(LMAX)
ENDIF
103 CONTINUE
DO 104 L1 = 1,3
L2 = LMAX + L1 -2
XOLD(L1) = XNEW(L2)
FOLD(L1) = FNEW(L2)
104 CONTINUE
IF (DX.LT.XACC)GOTO 10
100 CONTINUE
10 MAXTRI = XOLD(2)
RETURN
END
DOUBLE PRECISION FUNCTION RTBIS(FUNC, VAL,X1,X2,XACC)
DOUBLE PRECISION X1,X2,FUNC,VAL,XACC
DOUBLE PRECISION FMID,F,XMID,DX
INTEGER JMAX,J
PARAMETER (JMAX=100)
FMID=FUNC(X2)-VAL
F=FUNC(X1)-VAL
IF (F.LT.0.D0) THEN
RTBIS=X1
DX = X2-X1
ELSE
RTBIS=X2
DX = X1-X2
ENDIF
DO 11 J=1,JMAX
DX= 0.5D0*DX
XMID=RTBIS + DX
FMID=FUNC(XMID)-VAL
IF (FMID.LE.0.D0)RTBIS = XMID
IF (DABS(DX).LT.XACC) RETURN
11 CONTINUE
END
DOUBLE PRECISION FUNCTION RTSEC(SFN1,S,TLOW,TUP,XACC)
DOUBLE PRECISION TLOW,TUP,XACC,XL,SFN1,S
DOUBLE PRECISION FL,F,SWAP,DX
INTEGER MAXIT,H
PARAMETER (MAXIT=30)
FL=SFN1(TLOW)-S
F=SFN1(TUP)-S
IF (ABS(FL).LT.ABS(F)) THEN
RTSEC=TLOW
XL=TUP
SWAP=FL
FL=F
F=SWAP
ELSE
XL=TLOW
RTSEC=TUP
ENDIF
DO 11 H=1,MAXIT
DX=(XL-RTSEC)*F/(F-FL)
XL=RTSEC
FL=F
RTSEC=RTSEC+DX
F=SFN1(RTSEC) -S
IF (DABS(DX).LT.XACC) RETURN
11 CONTINUE
END
DOUBLE PRECISION FUNCTION RD(X,Y,Z,ERRTOL,IERR)
INTEGER IERR
DOUBLE PRECISION C1,C2,C3,C4,EA,EB,EC,ED,EF,EPSLON,ERRTOL,LAMDA
DOUBLE PRECISION LOLIM,MU,POWER4,SIGMA,S1,S2,UPLIM,X,XN,XNDEV
DOUBLE PRECISION XNROOT,Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
DATA LOLIM/1.0D-25/,UPLIM/4.5D+21/
IF (DMIN1(X,Y).LT.0.D0) GOTO 100
IF (DMIN1(X+Y,Z).LT.LOLIM) GOTO 100
IF (DMAX1(X,Y,Z).LE.UPLIM) GOTO 112
100 PRINT *,'INVALID ARGUMENTS PASSED TO RD'
IERR=1
GOTO 124
112 IERR=0
XN=X
YN=Y
ZN=Z
SIGMA=0.D0
POWER4=1.D0
116 MU=(XN+YN+3.D0*ZN)*0.2D0
XNDEV=(MU-XN)/MU
YNDEV=(MU-YN)/MU
ZNDEV=(MU-ZN)/MU
EPSLON=DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
IF (EPSLON.LT.ERRTOL) GOTO 120
XNROOT=DSQRT(XN)
YNROOT=DSQRT(YN)
ZNROOT=DSQRT(ZN)
LAMDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
SIGMA=SIGMA+POWER4/(ZNROOT*(ZN+LAMDA))
POWER4=POWER4*0.25D0
XN=(XN+LAMDA)*0.25D0
YN=(YN+LAMDA)*0.25D0
ZN=(ZN+LAMDA)*0.25D0
GOTO 116
120 C1=3.D0/14.D0
C2=1.D0/6.D0
C3=9.D0/22.D0
C4=3.D0/26.D0
EA=XNDEV*YNDEV
EB=ZNDEV*ZNDEV
EC=EA-EB
ED=EA-6.D0*EB
EF=ED+EC+EC
S1=ED*(-C1+0.25D0*C3*ED-1.5D0*C4*ZNDEV*EF)
S2=ZNDEV*(C2*EF+ZNDEV*(-C3*EC+ZNDEV*C4*EA))
RD=3.D0*SIGMA+POWER4*(1.D0+S1+S2)/(MU*DSQRT(MU))
124 RETURN
END
DOUBLE PRECISION FUNCTION RF(X,Y,Z,ERRTOL,IERR)
INTEGER IERR
DOUBLE PRECISION C1,C2,C3,E2,E3,EPSLON,ERRTOL,LAMBDA
DOUBLE PRECISION LOLIM,MU,S,UPLIM,X,XN,XNDEV,XNROOT
DOUBLE PRECISION Y,YN,YNDEV,YNROOT,Z,ZN,ZNDEV,ZNROOT
DATA LOLIM/1.5D-38/, UPLIM/ 3.0D+37/
IF (DMIN1(X,Y,Z).LT.0.D0) GOTO 100
IF (DMIN1(X+Y,X+Z,Y+Z).LT.LOLIM) GOTO 100
IF (DMAX1(X,Y,Z).LE.UPLIM) GOTO 112
100 PRINT *,'INVALID ARGUMENTS PASSED TO RF'
IERR=1
GOTO 124
112 IERR=0
XN=X
YN=Y
ZN=Z
116 MU=(XN+YN+ZN)/3.D0
XNDEV=2.D0-(MU+XN)/MU
YNDEV=2.D0-(MU+YN)/MU
ZNDEV=2.D0-(MU+ZN)/MU
EPSLON=DMAX1(DABS(XNDEV),DABS(YNDEV),DABS(ZNDEV))
IF (EPSLON.LT.ERRTOL) GOTO 120
XNROOT=DSQRT(XN)
YNROOT=DSQRT(YN)
ZNROOT=DSQRT(ZN)
LAMBDA=XNROOT*(YNROOT+ZNROOT)+YNROOT*ZNROOT
XN=(XN+LAMBDA)*0.25D0
YN=(YN+LAMBDA)*0.25D0
ZN=(ZN+LAMBDA)*0.25D0
GOTO 116
120 C1=1.D0/24.D0
C2=3.D0/44.D0
C3=1.D0/14.D0
E2=XNDEV*YNDEV-ZNDEV*ZNDEV
E3=XNDEV*YNDEV*ZNDEV
S=1.D0+(C1*E2-0.1D0-C2*E3)*E2+C3*E3
RF=S/DSQRT(MU)
124 RETURN
END
DOUBLE PRECISION FUNCTION DACOSH(X)
DOUBLE PRECISION X
DACOSH = DLOG (X + DSQRT(X*X-1.D0))
RETURN
END
Listing 1c - ISHP4
PROGRAM ISHP4
C ISHP version 4.06 - AUTOMATIC RUNNER
C
C Written by Paul F. Johnson, Feb. 1999 for his project
C
C This version has a few bug fixes and a small change to the
C save routine (adds the .csv extension)
C
C Calculation and comparison of the difference between exact
C scattering, ISHP, and the Porod's law approximation, ISHPLOC
C
C The main calculation, Porod and Exact routines have all been
C changed to make debugging a lot simpler.
C
C The results are stored as a CSV file.
C
C This version has been upgraded from ISHP3.0x to make debugging
C a hell of a lot simpler to perform and also to make the
C program more flexible for the future.
C
C It outputs the four basic types of graph :-
C
C 1. I(Q) vs Q (basic raw plot)
C 2. Ln I(Q) vs Q^2 (Guinier regime - duff for approx forms)
C 3. Log I(Q) vs Q (really good for low intensities of I(Q))
C 4. Log I(Q) vs Log I(Q) (power law / fractal scattering)
C
C Here's hoping!
COMMON/ISHP/QVEC,NQ,R0,MI,LI,GP,RPC
COMMON/FILER/FNAME
COMMON/CHANGE/XVEC,YVEC
REAL QVEC,XVEC,YVEC
CHARACTER*12 FNAME
CHARACTER*2 YN
DIMENSION QVEC(1000),XVEC(1000),YVEC(1000)
REAL R0,DQ,QMAX
INTEGER NQ,MI,LI,GP,RPC
INTEGER LOOP
1 WRITE (*,*)
WRITE (*,'('' PROGRAM ISHP4.06 BY PAUL F. JOHNSON, FEB 1999'')')
WRITE (*,*)
WRITE (*,'('' THIS PROGRAM ALLOWS THE INVESTIGATION OF SMALL '')')
WRITE (*,'('' ANGLE SCATTERING OF I(Q) FOR THE EXACT OR POROD'')')
WRITE (*,'('' APPROXIMATION.'')')
WRITE (*,*)
WRITE (*,'('' THE Q TERM CAN BE EITHER THE APPROXIMATE FOR OR'')')
WRITE (*,'('' THE Q**4 FORM. THE FINAL FIGURE. THE OUTPUT '')')
WRITE (*,'('' FILES CAN ALSO BE OF ONE OF FOUR FORMS (ALL'')')
WRITE (*,'('' FILES SAVED AS CSV)'')')
WRITE (*,*)
WRITE (*,'('' THIS VERSION WILL AUTOMATICALLY PERFORM ALL THE'')')
WRITE (*,'('' CALCULATIONS, JUST ASKING YOU FOR THE FILENAMES'')')
WRITE (*,*)
WRITE (*,'('' THE FOUR FORMS ARE : '')')
WRITE (*,'('' 1. I(Q) vs Q (THE NORMAL PLOT)'')')
WRITE (*,'('' 2. LN I(Q) vs Q**2 - THE GUINIER REGIME'')')
WRITE (*,'('' 3. LOG I(Q) vs Q - VERY GOOD FOR LOW INTENSITY'')')
WRITE (*,'('' I(Q)'')')
WRITE (*,'('' 4. LOG I(Q) vs LOG I(Q) - FOR POWER LAW / '')')
WRITE (*,'('' FRACTAL SCATTERING'')')
WRITE (*,*)
WRITE (*,'('' PLEASE ENTER QMAX '',$)')
READ *,QMAX
2 WRITE (*,'('' PLEASE ENTER THE NUMBER OF POINTS '',$)')
READ *,NQ
IF (NQ.LT.2) THEN
WRITE (*,'('' YOU NEED MORE THAN 2 POINTS.'')')
GOTO 2
ENDIF
IF (NQ.GT.1000) THEN
WRITE (*,'('' NO MORE THAN 1000 POINTS.'')')
GOTO 1
ENDIF
C Next, we need r0, the radius
4 WRITE (*,'('' PLEASE ENTER THE RADIUS OF THE SPHERE '',$)')
READ *,R0
IF (R0.LE.0) THEN
WRITE (*,'('' YOU NEED A RADIUS GREATER THAN ZERO!'')')
GOTO 4
ENDIF
WRITE (*,*)
WRITE (*,'('' ARE YOU USING AN ACORN MACHINE (Y/N)? '',$)')
READ (*,'(A)') YN
IF (YN.EQ.'Y'.OR.YN.EQ.'y') THEN
RPC=1
ELSE
RPC=2
ENDIF
C Okay, we'll now set QVEC, once and once only. Okay Trev :-)
DQ=QMAX/NQ
DO 5 LOOP=1,NQ
QVEC(LOOP)=((LOOP-1)+0.5E0)*DQ
5 CONTINUE
C Now that's out the way, let's do the loops.
C Variables : MI - Plot counter
C LI - Q / Q**4 flag
C GP - Guinier or Porod flag (1 = Guinier)
C This loop should be enough
LI=1
GP=1
WRITE (*,*)
6 DO 10 MI=1,4
IF (GP.EQ.1) THEN
IF (LI.EQ.1) THEN
WRITE(*,'(''EXACT FORM, EQUATION '',$)')
ELSE
WRITE(*,'(''EXACT FORM, Q**4, EQUATION '',$)')
ENDIF
ENDIF
C ELSE
IF (GP.EQ.2) THEN
IF (LI.EQ.1) THEN
WRITE(*,'(''POROD FORM, EQUATION '',$)')
ELSE
WRITE(*,'(''POROD FORM, Q**4, EQUATION '',$)')
ENDIF
ENDIF
WRITE(*,*) MI
CALL CALC()
CALL FILR()
CALL SAVER
10 CONTINUE
IF (LI.EQ.1.AND.GP.EQ.1) THEN
LI=2
GOTO 6
ENDIF
IF (LI.EQ.2.AND.GP.EQ.1) THEN
GP=GP+1
LI=1
GOTO 6
ENDIF
IF (LI.EQ.1.AND.GP.EQ.2) THEN
LI=2
GOTO 6
ENDIF
11 WRITE (*,*)
WRITE (*,'(''DO YOU WANT ANOTHER RUN (Y/N) ?'',$)')
READ (*,'(A)') YN
IF (YN.EQ.'Y'.OR.YN.EQ.'y') GOTO 1
END
C SUBROUTINE CALC.
C
C This routine gets rid of both the old subroutines DOEXACT,
C DOPOROD, EXACT() and POROD() in one easy to digest lump.
C The code is pretty straight forward.
SUBROUTINE CALC()
COMMON/ISHP/QVEC,NQ,R0,MI,LI,GP,RPC
COMMON/FILER/FNAME
COMMON/CHANGE/XVEC,YVEC
REAL QVEC,XVEC,YVEC,TEMP
CHARACTER*12 FNAME
DIMENSION QVEC(1000),XVEC(1000),YVEC(1000)
REAL R0,DQ,QMAX
INTEGER NQ,MI,LI,GP,RPC
INTEGER LOOP
REAL EISHP,PISHP
DO 10 LOOP=1,NQ
C Okay, let's deal with the X axis - XVEC
IF (MI.EQ.1.OR.MI.EQ.3) XVEC(LOOP)=QVEC(LOOP)
IF (MI.EQ.2) XVEC(LOOP)=QVEC(LOOP)**2
IF (MI.EQ.4) XVEC(LOOP)=LOG(QVEC(LOOP))
C Okay, let's have a look at the YVEC
IF (GP.EQ.1) THEN
TEMP=EISHP(QVEC(LOOP),R0,LI)
IF (MI.EQ.1) YVEC(LOOP)=TEMP
IF (MI.EQ.2) YVEC(LOOP)=ALOG(TEMP)
IF (MI.EQ.3.OR.MI.EQ.4) YVEC(LOOP)=LOG(TEMP)
ELSE
TEMP=PISHP(QVEC(LOOP),R0,LI)
IF (MI.EQ.1) YVEC(LOOP)=TEMP
IF (MI.EQ.2) YVEC(LOOP)=ALOG(TEMP)
IF (MI.EQ.3.OR.MI.EQ.4) YVEC(LOOP)=LOG(TEMP)
ENDIF
10 CONTINUE
RETURN
END
REAL FUNCTION EISHP(Q1,R0,LI)
C This function calculates the exact scattering Ishp for a
C given Q and R0.
REAL QR0,FSPH,PI,Q1,R0,PASS
INTEGER LI
PI=3.1415927E0
QR0=Q1*R0
FSPH=4*PI*(Q1**(-3))*(SIN(QR0)-((QR0)*COS(QR0)))
EISHP=FSPH*FSPH
IF (LI.EQ.2) THEN
PASS=EISHP*(Q1**4)
EISHP=PASS
ENDIF
RETURN
END
REAL FUNCTION PISHP(Q1,R0,LI)
C This function calculates the scattering accoring to Porod
REAL QR0,ASPH,PI,Q1,R0,PASS
INTEGER LI
PI=3.1415927E0
QR0=Q1*R0
ASPH=4*PI*R0*R0
PISHP=2*PI*ASPH*(Q1**(-4))*(1+QR0**(-2))
IF (LI.EQ.2) THEN
PASS=PISHP*(Q1**4)
PISHP=PASS
ENDIF
RETURN
END
SUBROUTINE FILR()
C This subroutine assigns the filename and adds an extension
C depending on the platform (i.e. PCs get .CSV added)
COMMON/ISHP/QVEC,NQ,R0,MI,LI,GP,RPC
COMMON/FILER/FNAME
COMMON/CHANGE/XVEC,YVEC
REAL QVEC,XVEC,YVEC
CHARACTER*12 FNAME
CHARACTER*8 NAME
CHARACTER*4 EXT
DIMENSION QVEC(1000),XVEC(1000),YVEC(1000)
REAL R0,DQ,QMAX
INTEGER NQ,MI,LI,GP,RPC
EXT='.CSV'
WRITE (*,*)
IF (RPC.EQ.1) THEN
WRITE (*,'('' FILENAME (10 CHARS MAX.) '',$)')
READ (*,'(A)') FNAME
ELSE
WRITE (*,'('' FILENAME (8 CHARS MAX.) '',$)')
READ (*,'(A)') NAME
ENDIF
IF (RPC.EQ.2) Fname=NAME//EXT
RETURN
END
SUBROUTINE SAVER()
COMMON/ISHP/QVEC,NQ,R0,MI,LI,GP,RPC
COMMON/FILER/FNAME
COMMON/CHANGE/XVEC,YVEC
REAL QVEC,XVEC,YVEC
CHARACTER*12 FNAME
DIMENSION QVEC(1000),XVEC(1000),YVEC(1000)
REAL R0,DQ,QMAX
INTEGER NQ,MI,LI,GP,RPC,L
1 OPEN (4,FILE=FNAME,FORM='FORMATTED',ERR=999)
WRITE (4,*) 'XVEC(L)',',','YVEC(L)'
DO 2 L=1,NQ
WRITE (4,*) XVEC(L),',',YVEC(L)
2 CONTINUE
CLOSE (4)
RETURN
999 WRITE (*,'('' FILER ERROR.'')')
GOTO 1
END