Module BARKLEM PRIVATE PUBLIC :: RETCROSS CONTAINS subroutine RETCROSS(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) !********************************************************************** ! Returns the cross section at 10000 m/s from tables for neutrals etc ! returns IFAIL=1 if not on the tables or some other error, else 0 ! Paul Barklem Dec 1997 !********************************************************************** ! ! NSupp,NSlow = upper and lower state effective principal quantum numbers ! Lupp,Llow = angular momentum quantum numbers ! CROSS = hydrogen broadening cross-section for v=10000m/s ! ALPHA = velocity parameter ! implicit none real*4 NSupp,NSlow,CROSS,ALPHA integer Lupp,Llow,IFAIL,CHECK,TABLE ! CHECK=ABS(Lupp-Llow) if (CHECK.ne.1) then IFAIL=1 print *,'#######Forbidden transition!' return end if ! TABLE=Llow+Lupp if (TABLE.eq.1) & call SP(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) if (TABLE.eq.3) & call PD(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) if (TABLE.eq.5) & call DF(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) if ((TABLE.ne.1).and.(TABLE.ne.3).and.(TABLE.ne.5)) & then IFAIL=1 print *,'#######Transition not covered by current data' return end if ! return end subroutine ! subroutine SP(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) !********************************************************************** ! Returns the cross section at 10000 m/s from s-p table ! returns IFAIL=1 if not on the tables, else 0 ! Paul Barklem Dec 1997 !********************************************************************** ! implicit none real*4 NSupp,NSlow,CROSS,ALPHA,NSS,NSP, & CS(21,18),AL(21,18),CS2(21,18),AL2(21,18), & NSSG(21),NSPG(18) integer Lupp,Llow,IFAIL,I,J character*70 COMMENTS ! if (Lupp.eq.1) then NSP=NSupp NSS=NSlow else NSP=NSlow NSS=NSupp end if ! if ((NSS.gt.3.).or.(NSS.lt.1.).or.(NSP.gt.3.).or.(NSP.lt.1.3)) & then IFAIL=1 print *,'#######Outside range of tables! (s-p)' return end if ! ! read in the table data ! 3 lines of comments preceed ! open(15,file='data/barklem/spdata.dat',form='formatted', & status='old') do I=1,3 read(15,'(A70)') COMMENTS ! print *,COMMENTS end do ! do I=1,21 read(15,*) (CS(I,J),J=1,18) end do read(15,'(A70)') COMMENTS ! a single comment line between data do I=1,21 read(15,*) (AL(I,J),J=1,18) end do ! do I=1,21 NSSG(I)=1.0+(I-1)*0.1 end do do I=1,18 NSPG(I)=1.3+(I-1)*0.1 end do ! ! setup second derivative table for spline ! call SPLIE2(NSSG,NSPG,CS,21,18,CS2) call SPLIE2(NSSG,NSPG,AL,21,18,AL2) ! ! run bicubic spline interpolation ! call SPLIN2(NSSG,NSPG,CS,CS2,21,18,NSS,NSP,CROSS) call SPLIN2(NSSG,NSPG,AL,AL2,21,18,NSS,NSP,ALPHA) ! close(15) return end subroutine ! subroutine PD(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) !********************************************************************** ! Returns the cross section at 10000 m/s from p-d table ! returns IFAIL=1 if not on the tables, else 0 ! Paul Barklem Dec 1997 !********************************************************************** ! implicit none real*4 NSupp,NSlow,CROSS,ALPHA,NSD,NSP, & CS(18,18),AL(18,18),CS2(18,18),AL2(18,18), & NSPG(18),NSDG(18) integer Lupp,Llow,IFAIL,I,J character*70 COMMENTS ! if (Lupp.eq.2) then NSD=NSupp NSP=NSlow else NSD=NSlow NSP=NSupp end if ! if ((NSP.gt.3.).or.(NSP.lt.1.3).or.(NSD.gt.4.).or.(NSD.lt.2.3)) & then IFAIL=1 print *,'#######Outside range of tables! (p-d)' return end if ! ! read in the table data ! 3 lines of comments preceed ! open(15,file='data/barklem/pddata.dat',form='formatted', & status='old') do I=1,3 read(15,'(A70)') COMMENTS ! print *,COMMENTS end do ! do I=1,18 read(15,*) (CS(I,J),J=1,18) end do read(15,'(A70)')COMMENTS ! a single comment line between data do I=1,18 read(15,*) (AL(I,J),J=1,18) end do ! ! make arrays of the nstar values ! do I=1,18 NSPG(I)=1.3+(I-1)*0.1 NSDG(I)=2.3+(I-1)*0.1 end do ! ! setup second derivative table for spline ! call SPLIE2(NSPG,NSDG,CS,18,18,CS2) call SPLIE2(NSPG,NSDG,AL,18,18,AL2) ! ! run bicubic spline interpolation ! call SPLIN2(NSPG,NSDG,CS,CS2,18,18,NSP,NSD,CROSS) call SPLIN2(NSPG,NSDG,AL,AL2,18,18,NSP,NSD,ALPHA) ! ! close(15) return end subroutine ! subroutine DF(NSlow,NSupp,Llow,Lupp,CROSS,ALPHA,IFAIL) !********************************************************************** ! Returns the cross section at 10000 m/s from d-f table ! returns IFAIL=1 if not on the tables, else 0 ! Paul Barklem Dec 1997 !********************************************************************** ! implicit none real*4 NSupp,NSlow,CROSS,ALPHA,NSD,NSF, & CS(18,18),AL(18,18),CS2(18,18),AL2(18,18), & NSDG(18),NSFG(18) integer Lupp,Llow,IFAIL,I,J character*70 COMMENTS ! if (Lupp.eq.3) then NSF=NSupp NSD=NSlow else NSF=NSlow NSD=NSupp end if ! if ((NSD.gt.4.).or.(NSD.lt.2.3).or.(NSF.gt.5.).or.(NSF.lt.3.3)) & then IFAIL=1 print *,'#######Outside range of tables! (d-f)' return end if ! ! read in the table data ! 3 lines of comments preceed ! open(15,file='data/barklem/dfdata.dat',form='formatted', & status='old') do I=1,3 read(15,'(A70)') COMMENTS ! print *,COMMENTS end do ! do I=1,18 read(15,*) (CS(I,J),J=1,18) ! print *,(CS(I,J),J=1,18) end do read(15,'(A70)')COMMENTS ! a single comment line between data do I=1,18 read(15,*) (AL(I,J),J=1,18) ! print *,(AL(I,J),J=1,18) end do ! ! make arrays of the nstar values ! do I=1,18 NSDG(I)=2.3+(I-1)*0.1 NSFG(I)=3.3+(I-1)*0.1 end do ! ! setup second derivative table for spline ! call SPLIE2(NSDG,NSFG,CS,18,18,CS2) call SPLIE2(NSDG,NSFG,AL,18,18,AL2) ! ! run bicubic spline interpolation ! call SPLIN2(NSDG,NSFG,CS,CS2,18,18,NSD,NSF,CROSS) call SPLIN2(NSDG,NSFG,AL,AL2,18,18,NSD,NSF,ALPHA) ! ! close(15) return end subroutine ! subroutine WIDTH(CROSS,ALPHA,TEMP,A,HALFWIDTH) !********************************************************************** ! Computes the linewidth (half width at half maxima) given ! the cross-section and velocity parameter and tempertature ! Paul Barklem Dec 1997 !********************************************************************** ! ! CROSS = cross section in atomic units at 10000m/s ! ALPHA = velocity parameter ! TEMP = temperature in Kelvin ! HALFWIDTH = half width per unit hydrogen atom density (rad /s cm^-3) ! A = Atomic Mass of the absorbing atom eg. A(Mg)=24.32 ! RMASS = reduced mass ! MEANVEL = mean kinetic particle velocity in m/s ! implicit none real*4 CROSS,ALPHA,TEMP,HALFWIDTH,MEANVEL,RMASS,CROSSMEAN, & CROSSM,A,PI,K,M0,A0 ! real*8 DGAMMA ! parameter (PI=3.141592653) parameter (K=1.380658E-23) !boltzmanns constant J/K parameter (M0=1.660540E-27) !unit atomic mass kg (Carbon 12) parameter (A0=5.29177249E-11) !bohr radius m ! CROSSM=CROSS*A0*A0 !cross section to m^2 RMASS=M0/(1./1.008+1./A) !calculate reduced mass MEANVEL=SQRT(8.0*K*TEMP/PI/RMASS) !array of mean velocities ! ! work out cross section at mean perturber velocity for this Temp ! CROSSMEAN= CROSSM*((MEANVEL/10000.0)**(-ALPHA)) ! ! the half-half width per unit perturber density m^3 rad/s ! HALFWIDTH=(4.0/PI)**(ALPHA/2.0) * & SNGL(DGAMMA((4.0-ALPHA)/2.d0) )* & MEANVEL*CROSSMEAN ! HALFWIDTH=HALFWIDTH * 1.e6 ! to rad/s cm^3 ! return end subroutine ! !********************************************************************** ! The next routine is from the NETLIB archive ! Computes the Gamma function !********************************************************************** DOUBLE PRECISION FUNCTION DGAMMA(X) !---------------------------------------------------------------------- ! From NETLIB, group SPECFUN ! 'AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL FUNCTIONS', ! LECTURE NOTES IN MATHEMATICS, 506, NUMERICAL ANALYSIS DUNDEE, ! 1975, G. A. WATSON (ED.), SPRINGER VERLAG, BERLIN, 1976. !******************************************************************* ! ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS ! ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0 + EPS .GT. 1.0 ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION ! GAMMA(XBIG) = XINF. ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1/XMININ IS MACHINE REPRESENTABLE. ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER. ! ! ! ERROR RETURNS ! ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED ! TO BE FREE OF UNDERFLOW AND OVERFLOW. ! ! OTHER SUBPROGRAMS REQUIRED (DOUBLE PRECISION VERSION) ! ! DBLE,DEXP,DLOG,DSIN,FLOAT,IFIX,SNGL ! double precision C,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI, & SUM,TWELVE,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO integer I,J,N logical PARITY dimension C(7),P(8),Q(8) !---------------------------------------------------------------------- ! MATHEMATICAL CONSTANTS !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,ZERO/1.0D0,0.5D0,12.0D0,0.0D0/ DATA SQRTPI/0.9189385332046727417803297D0/ DATA PI/3.1415926535897932384626434D0/ !---------------------------------------------------------------------- ! MACHINE DEPENDENT PARAMETERS !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/34.844D0,5.883D-39,2.776D-17/ DATA XINF/1.7014D38/ !---------------------------------------------------------------------- ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX ! APPROXIMATION OVER (1,2). !---------------------------------------------------------------------- DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, & -3.79804256470945635097577D+2,6.29331155312818442661052D+2, & 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, & -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, & -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, & 2.25381184209801510330112D+4,4.75584627752788110767815D+3, & -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- DATA C/-1.910444077728D-03,8.4171387781295D-04, & -5.952379913043012D-04,7.93650793500350248D-04, & -2.777777777777681622553D-03,8.333333333333333331554247D-02, & 5.7083835261D-03/ !---------------------------------------------------------------------- PARITY = .FALSE. FACT = ONE N = 0 Y = X IF (Y .GT. ZERO) GO TO 200 !---------------------------------------------------------------------- ! ARGUMENT IS NEGATIVE !---------------------------------------------------------------------- Y = -X J = IFIX(SNGL(Y)) RES = Y - DBLE(FLOAT(J)) IF (RES .EQ. ZERO) GO TO 700 IF (J .NE. (J/2)*2) PARITY = .TRUE. FACT = -PI / DSIN(PI*RES) Y = Y + ONE !---------------------------------------------------------------------- ! ARGUMENT IS POSITIVE !---------------------------------------------------------------------- 200 IF (Y .LT. EPS) GO TO 650 IF (Y .GE. TWELVE) GO TO 300 Y1 = Y IF (Y .GE. ONE) GO TO 210 !---------------------------------------------------------------------- ! 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- Z = Y Y = Y + ONE GO TO 250 !---------------------------------------------------------------------- ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY !---------------------------------------------------------------------- 210 N = IFIX(SNGL(Y)) - 1 Y = Y - DBLE(FLOAT(N)) Z = Y - ONE !---------------------------------------------------------------------- ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 !---------------------------------------------------------------------- 250 XNUM = ZERO XDEN = ONE DO 260 I = 1, 8 XNUM = (XNUM + P(I)) * Z XDEN = XDEN * Z + Q(I) 260 CONTINUE RES = XNUM / XDEN + ONE IF (Y .EQ. Y1) GO TO 900 IF (Y1 .GT. Y) GO TO 280 !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- RES = RES / Y1 GO TO 900 !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 !---------------------------------------------------------------------- 280 DO 290 I = 1, N RES = RES * Y Y = Y + ONE 290 CONTINUE GO TO 900 !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- 300 IF (Y .GT. XBIG) GO TO 700 YSQ = Y * Y SUM = C(7) DO 350 I = 1, 6 SUM = SUM / YSQ + C(I) 350 CONTINUE SUM = SUM/Y - Y + SQRTPI SUM = SUM + (Y-HALF)*DLOG(Y) RES = DEXP(SUM) GO TO 900 !---------------------------------------------------------------------- ! ARGUMENT .LT. EPS !---------------------------------------------------------------------- 650 IF (Y .LT. XMININ) GO TO 700 RES = ONE / Y GO TO 900 !---------------------------------------------------------------------- ! RETURN FOR SINGULARITIES, EXTREME ARGUMENTS, ETC. !---------------------------------------------------------------------- 700 DGAMMA = XINF print *,'function DGAMMA has extreme argument :' print *,'argument is : ',X print *,'function returns machine upper bound ' GO TO 950 !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- 900 IF (PARITY) RES = -RES IF (FACT .NE. ONE) RES = FACT / RES DGAMMA = RES 950 RETURN ! ---------- LAST CARD OF GAMMA ---------- END FUNCTION ! !********************************************************************** ! The following are Numerical Recipes routines ! W. Press, S. Teukolsky, W. Vetterling and B. Flannery ! " Numerical Recipes in Fortran: The Art of Scientific ! Computing" Second Edition, Cambridge Univ. Press !********************************************************************** ! SUBROUTINE SPLIE2(X1A,X2A,YA,M,N,Y2A) PARAMETER (NN=100) DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN) DO 13 J=1,M DO 11 K=1,N YTMP(K)=YA(J,K) 11 CONTINUE CALL SPLINE(X2A,YTMP,N,1.E30,1.E30,Y2TMP) DO 12 K=1,N Y2A(J,K)=Y2TMP(K) 12 CONTINUE 13 CONTINUE RETURN END SUBROUTINE SUBROUTINE SPLIN2(X1A,X2A,YA,Y2A,M,N,X1,X2,Y) PARAMETER (NN=100) DIMENSION X1A(M),X2A(N),YA(M,N),Y2A(M,N),YTMP(NN),Y2TMP(NN), & YYTMP(NN) DO 12 J=1,M DO 11 K=1,N YTMP(K)=YA(J,K) Y2TMP(K)=Y2A(J,K) 11 CONTINUE CALL SPLINT(X2A,YTMP,Y2TMP,N,X2,YYTMP(J)) 12 CONTINUE CALL SPLINE(X1A,YYTMP,M,1.E30,1.E30,Y2TMP) CALL SPLINT(X1A,YYTMP,Y2TMP,M,X1,Y) RETURN END SUBROUTINE SUBROUTINE SPLINE(X,Y,N,YP1,YPN,Y2) PARAMETER (NMAX=100) DIMENSION X(N),Y(N),Y2(N),U(NMAX) IF (YP1.GT..99E30) THEN Y2(1)=0. U(1)=0. ELSE Y2(1)=-0.5 U(1)=(3./(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) ENDIF DO 11 I=2,N-1 SIG=(X(I)-X(I-1))/(X(I+1)-X(I-1)) P=SIG*Y2(I-1)+2. Y2(I)=(SIG-1.)/P U(I)=(6.*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) & /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 11 CONTINUE IF (YPN.GT..99E30) THEN QN=0. UN=0. ELSE QN=0.5 UN=(3./(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) ENDIF Y2(N)=(UN-QN*U(N-1))/(QN*Y2(N-1)+1.) DO 12 K=N-1,1,-1 Y2(K)=Y2(K)*Y2(K+1)+U(K) 12 CONTINUE RETURN END SUBROUTINE SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y) DIMENSION XA(N),YA(N),Y2A(N) KLO=1 KHI=N 1 IF (KHI-KLO.GT.1) THEN K=(KHI+KLO)/2 IF(XA(K).GT.X)THEN KHI=K ELSE KLO=K ENDIF GOTO 1 ENDIF H=XA(KHI)-XA(KLO) ! TM remplace PAUSE par STOP IF (H.EQ.0.) STOP 'Bad XA input.' A=(XA(KHI)-X)/H B=(X-XA(KLO))/H Y=A*YA(KLO)+B*YA(KHI)+ & ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6. RETURN END SUBROUTINE !---------------------------------------------------------------------- End Module BARKLEM !********************************************************************** !********************************************************************** MODULE BURGESS_TULLY ! (A&A Burgess & Tully 1992) ! Calcul des moyennes thermiques des forces de collisions ! pour des ions positifs ! Ex : CaII, MgII ... ! ! Upsil = f(Spline, T, transition) ! Implicit None ! PRIVATE PUBLIC :: UPSIL, OMEGA ! CONTAINS ! REAL FUNCTION SPLINE(P1,P2,P3,P4,P5,X) Real, Intent(In) :: P1, P2, P3, P4, P5, X ! 5-point spline interpolation of Y(X), for X in the range (0,1). ! Input: ! Knot values Pl=Y(O), P2=Y(1/4), P3=Y(1/2), P4=Y(3/4), P5=Y(1} ! and X. ! Output: C SPLINE=Y(X) . ! code: Real :: S, S2, S3, S4, X0, T0, T1, T2, T3 S=1.0/30.0 S2=32.0*S*(19.0*P1-43.0*P2+30.0*P3-7.0*P4+P5) S3=160.0*S*(-P1+7.0*P2-12.0*P3+7.0*P4-P5) S4=32.0*S*(P1-7.0*P2+30.0*P3-43.0*P4+19.0*P5) IF(X.GT.0.25) GO TO 1 X0=X-0.125 T3=0.0 T2=0.5*S2 T1=4.0*(P2-P1) T0=0.5*(P1+P2)-0.015625*T2 GO TO 4 ! 1 IF(X.GT.0.5) GO TO 2 X0=X-0.375 T3=20.0*S*(S3-S2) T2=0.25*(S2+S3) T1=4.0*(P3-P2)-0.015625*T3 T0=0.5*(P2+P3)-0.015625*T2 GO TO 4 ! 2 IF(X.GT.0.75) GO TO 3 X0=X-0.625 T3=20.0*S*(S4-S3) T2=0.25*(S3+S4) T1=4.0*(P4-P3)-0.015625*T3 T0=0.5*(P3+P4)-0.015625*T2 GO TO 4 ! 3 X0=X-0.875 T3=0.0 T2=0.5*S4 T1=4.0*(P5-P4) T0=0.5*(P4+P5)-0.015625*T2 ! 4 SPLINE=T0+X0*(T1+X0*(T2+X0*T3)) END FUNCTION SPLINE !CC !CC !CC REAL FUNCTION UPSIL(K,EIJ,C,P1,P2,P3,P4,P5,T) Real, Intent(In) :: EIJ, C, P1, P2, P3, P4, P5, T Integer, Intent(In) :: K ! Maxwell averaged collision strength, Upsilon. ! Input: ! K.transition type, ! EIJ.transition energy (Ryd.), ! !.abscissa scale parameter, ! Pl,P2,P3,P4,P5-spline knot values, ! T=electron temperature (Kelvin). ! Output: ! UPSIL=Upsilon. ! Code: ! Real :: E, X, Y E=ABS(T/(1.57888E5*EIJ)) IF ((K.EQ.1).OR.(K.EQ.4)) X=LOG((E+C)/C)/LOG(E+C) IF ((K.EQ.2).OR.(K.EQ.3)) X=E/(E+C) Y=SPLINE(P1,P2,P3,P4,P5,X) IF (K.EQ.1) Y=Y*LOG(E+2.71828) IF (K.EQ.3) Y=Y/(E+1) IF (K.EQ.4) Y=Y*LOG(E+C) UPSIL=Y END FUNCTION UPSIL !CC !CC !CC REAL FUNCTION OMEGA(K,EIJ,C,P1,P2,P3,P4,P5,EJ) Real, intent(in) :: EIJ,C,P1,P2,P3,P4,P5,EJ Integer, intent(in) :: K ! Collision strength, Omega. !Input: ! K = transition type, ! EIJ= transitions energy (Ryd.), ! C = abscissa scale parameter, ! P1,P2,p3,P4,P5 = spline knot values, ! EJ = colliding electron energy after excitation (Ryd.). !Output: ! OMEGA= Omega. ! Code: ! Real :: E, X, Y E = ABS(EJ/EIJ) IF ((K==1) .OR. (K==4)) X = LOG((E+C)/C) / LOG(E+C) IF ((K.EQ.2).OR.(K.EQ.3)) X=E/(E+C) Y = SPLINE(P1,P2,P3,P4,P5,X) IF (K.EQ.1) Y=Y*LOG(E+2.71828) IF (K.EQ.3) Y=Y/((E+1)*(E+1)) IF (K.EQ.4) Y=Y*LOG(E+C) OMEGA = Y END FUNCTION OMEGA !CC END MODULE BURGESS_TULLY !********************************************************************** !********************************************************************** Module BESSEL ! ! Fonctions de Bessel modifiés de deuxième (ou 3ème) espèce d'ordre 0 et 1 ! ! Copyright 1996 Takuya OOURA ! PRIVATE PUBLIC :: dbesk0, dbesk1 CONTAINS ! ! Bessel K_0(x) function in double precision ! function dbesk0(x) implicit double precision (a - h, o - z) dimension a(0 : 15), b(0 : 111), c(0 : 134), d(0 : 39) data (a(i), i = 0, 15) / & 2.4307270476772195953d-12, 4.7091666363785304370d-10, & 6.7816861334344265568d-8, 6.7816840204737508252d-6, & 4.3402777777915334676d-4, 1.5624999999999872796d-2, & 2.5000000000000000448d-1, 9.9999999999999999997d-1, & 6.5878327432224993071d-12, 1.2083308769932888218d-9, & 1.6271062073716412046d-7, 1.4914719278555277887d-5, & 8.4603509071212245667d-4, 2.5248929932162333910d-2, & 2.7898287891460312491d-1, 1.1593151565841244874d-1 / data (b(i), i = 0, 13) / & -4.6430702971053162197d-13, 1.0377936059563728230d-11, & -1.0298475936392057807d-10, 5.3632747492333959219d-10, & -2.1674628861036068105d-10, -2.3316071545820437669d-8, & 2.2557819578691704059d-7, -9.2325694638587080009d-7, & -3.3569097781613661759d-6, 8.7355061305812582974d-5, & -6.8021202111645760475d-4, 2.7434654781323362319d-4, & 1.0031787169953909561d-1, 4.2102443824070833334d-1 / data (b(i), i = 14, 27) / & 4.1447451117883103686d-12, -3.4026589638576604315d-11, & 9.3398790624638977468d-12, 1.5184181750799852630d-9, & -1.1364911665083029464d-8, 2.0619457602095915719d-8, & 3.0431018037572243630d-7, -2.9749736264474555510d-6, & 8.0143661611467038568d-6, 8.0937525149549218398d-5, & -1.0356346549612699886d-3, 2.8534806627578638795d-3, & 9.7369634474060441807d-2, 3.2175066577856452683d-1 / data (b(i), i = 28, 41) / & 1.1170882570740727520d-13, -8.2865909408297066068d-11, & 9.4656678749191182763d-10, -3.5832019841847883380d-9, & -9.5017955656904252761d-9, 1.5200595674883329093d-7, & -3.8663262571356059980d-7, -3.3350340828235103499d-6, & 2.9359886663960844231d-5, -1.1266401822556801563d-5, & -1.2113572742435576205d-3, 6.3158973673701376253d-3, & 8.8291790250128171341d-2, 2.2833982383240512262d-1 / data (b(i), i = 42, 55) / & -3.2880638807053948433d-11, 4.3194884830465283512d-10, & -1.7455089683104033093d-9, -3.2437330799994764516d-9, & 4.7393655539139519778d-8, -1.1929265603456272466d-8, & -1.3177845881013419388d-6, 3.3873375636197969526d-6, & 3.2729835880668256625d-5, -1.8367283883002494561d-4, & -8.2830996454188084408d-4, 9.5512732229514251931d-3, & 7.2233832113719266702d-2, 1.4753187103603405298d-1 / data (b(i), i = 56, 69) / & 7.9998492614150860098d-11, -7.0257346702686139490d-10, & 7.8898821627084586270d-10, 1.1294796399671507085d-8, & -1.1360539648638059137d-8, -3.0346309115270564487d-7, & 3.2235585426189451721d-7, 8.3575612102298214948d-6, & -8.5169628089198208211d-6, -2.5740175232173357342d-4, & 1.2462734014689152770d-4, 1.0683232869192203450d-2, & 5.1515690033637395779d-2, 8.5465862953544883657d-2 / data (b(i), i = 70, 83) / & -8.6111506537356531608d-11, 5.1862926131024597823d-10, & 7.5884324949371110022d-10, -6.4011975813006767417d-9, & -4.1966181325111763156d-8, 9.1306285446881485314d-8, & 1.3573638315827954034d-6, 4.8683213252735694701d-7, & -3.8805424608710197066d-5, -1.1838986468688980610d-4,& 9.2796213947750964945d-4, 8.9611057737319027776d-3, & 3.1464453915862785606d-2, 4.4267648087536630780d-2 / data (b(i), i = 84, 97) / & 4.4400123834164610288d-11, -1.1411233140911074336d-10, & -8.8200670702467059830d-10, -1.9686735373323381456d-9, & 1.9921120728941773855d-8, 1.4543974418584834740d-7, & 1.8238418041265854754d-8, -4.5363700392899066037d-6, & -2.1688068222527688542d-5, 4.5496062166687034700d-5, & 1.0435238076080528284d-3, 5.8374528996419979931d-3, & 1.6611210710425455850d-2, 2.0756008367065750538d-2 / data (b(i), i = 98, 111) / & -6.5166519951106397214d-12, -5.8572182858788539580d-11, & 1.5550375065815375404d-10, 1.9526509484993563229d-9, & 9.2637123346818426594d-9, -1.4136471501812055943d-8, & -4.3024895710889717172d-7, -2.3235612243330592076d-6, & 4.0380616133862188804d-7, 9.2783767992909743602d-5, & 7.2964887597817095035d-4, 3.1316245282223273413d-3, & 7.8028233022066428316d-3, 9.0014807263791058095d-3 / data (c(i), i = 0, 14) / & 4.5161032649342790231d-11, -4.2774336988557091369d-11, & 6.0998467173896677777d-10, 1.9845167242599996944d-9, & 1.3097678767280215271d-8, 7.4505822268382641286d-8, & 4.2893920879106814989d-7, 2.3900851955655303104d-6, & 1.2533473009382380357d-5, 5.9693359063879871983d-5, & 2.4775070661087304580d-4, 8.5106703131389516508d-4, & 2.2500105115665788755d-3, 4.0446134454521634600d-3, & 3.6910983340425942762d-3 / data (c(i), i = 15, 29) / & 3.5732826433251464989d-12, -3.2906649482312266258d-12, & 7.0873811190464760555d-11, 2.9551320580484177120d-10, & 2.2776940472505079894d-9, 1.5175463612815010036d-8, & 9.9462487812170164133d-8, 6.1448757797853901100d-7, & 3.4869531882907360750d-6, 1.7615836644757657443d-5, & 7.6373536037879531886d-5, 2.7098571871205999668d-4, & 7.3399047381788927036d-4, 1.3439197177355085297d-3, & 1.2439943280131230863d-3 / data (c(i), i = 30, 44) / & 3.6343547836242523646d-13, 9.7997961751276137602d-14, & 1.0184692699811569047d-11, 6.1495184828957652064d-11, & 5.0238328349302602543d-10, 3.7498626376004337661d-9, & 2.6689445483857236307d-8, 1.7591899737346368084d-7, & 1.0486448307010701679d-6, 5.4986458466257148573d-6, & 2.4521456351751345323d-5, 8.8900942259143832228d-5, & 2.4483947714068300190d-4, 4.5418248688489693045d-4, & 4.2479574186923180694d-4 / data (c(i), i = 45, 59) / & 5.2460389348163395857d-14, 7.4802063026503503540d-14, & 2.0012201610651998417d-12, 1.4887306044735163359d-11, & 1.2946705414232940350d-10, 1.0391628915892803144d-9, & 7.8091180499677328456d-9, 5.3694223626907660084d-8, & 3.3063914804658509029d-7, 1.7776972424421486506d-6, & 8.0833148098458320202d-6, 2.9755556304448817780d-5, & 8.2945928349220642178d-5, 1.5536921180500112883d-4, & 1.4647070522281538711d-4 / data (c(i), i = 60, 74) / & 9.7531436733955514559d-15, 2.4084291220447154982d-14, & 4.7654956400897494468d-13, 4.0200949504810597783d-12, & 3.6726577109162191533d-11, 3.0939005665422637601d-10, & 2.4122848979784500179d-9, 1.7071884462645525505d-8, & 1.0752238955654933405d-7, 5.8844190041189462347d-7, & 2.7136083303224014597d-6, 1.0102477728604441135d-5, & 2.8420490721532571809d-5, 5.3637016379451944413d-5, & 5.0881312956459247572d-5 / data (c(i), i = 75, 89) / & 2.1732049868189377260d-15, 7.2720052142815590531d-15, & 1.2803083795536820100d-13, 1.1696825543787717167d-12, & 1.1083298191597132094d-11, 9.6536661252658773139d-11, & 7.7242553835198536397d-10, 5.5798366267110575620d-9, & 3.5721345296543414370d-8, 1.9806931547193682466d-7, & 9.2312964655319555313d-7, 3.4666258590861079959d-6, & 9.8224698307751177077d-6, 1.8648773453825584428d-5, & 1.7780062316167651812d-5 / data (c(i), i = 90, 104) / & 5.5012463763851934112d-16, 2.2254763392767319419d-15, & 3.7187669817701214965d-14, 3.5819585377733489628d-13, & 3.4866061263191556694d-12, 3.1101633450629652910d-11, & 2.5358235662235617663d-10, 1.8597629779492599046d-9, & 1.2052654739462999992d-8, 6.7501417351172136833d-8, & 3.1720052198654584574d-7, 1.1993651363602981832d-6, & 3.4179130317623363474d-6, 6.5208606745808860158d-6, & 6.2430205476536771454d-6 / data (c(i), i = 105, 119) / & 1.5225407517829491689d-16, 6.9834820025664405161d-16, & 1.1380182837138781431d-14, 1.1369488761077196511d-13, & 1.1291168681618466716d-12, 1.0250757630526871007d-11, & 8.4765287317253141514d-11, 6.2886627779402596211d-10, & 4.1142865598366029316d-9, 2.3223773435632014408d-8, & 1.0985095234166396934d-7, 4.1766260951820336228d-7, & 1.1958609263543792991d-6, 2.2907574647671878055d-6, & 2.2008253973114914005d-6 / data (c(i), i = 120, 134) / & 4.4863058691420695911d-17, 2.2437356594371819978d-16, & 3.6107964803015652759d-15, 3.7031193629853392081d-14, & 3.7341552790439784371d-13, 3.4355950129497564468d-12, & 2.8719942600171304499d-11, 2.1499646844509516453d-10, & 1.4171810843455227171d-9, 8.0501118772875784153d-9, & 3.8281889106330295876d-8, 1.4621673458431979989d-7, & 4.2029868696411098586d-7, 8.0785884122023473025d-7, & 7.7845438614204963209d-7 / data (d(i), i = 0, 7) / & -7.9737703860537066166d-14, 1.9543834380466766627d-12, & -4.7230794431646733538d-11, 1.4001773785771252004d-9, & -5.4864553020583098585d-8, 3.1601984250143742772d-6, & -3.3708783204090252161d-4, 1.6180215937964160437d-1 / data (d(i), i = 8, 15) / & -5.2593898374798632343d-14, 1.7725913926973236457d-12, & -4.6672234858122387294d-11, 1.3991653503828889207d-9, & -5.4863400156413929639d-8, 3.1601976099900075541d-6, & -3.3708783171335864627d-4, 1.6180215937958433760d-1 / data (d(i), i = 16, 23) / & -3.6135496189875398132d-14, 1.5466239429618130284d-12, & -4.5320259146602122624d-11, 1.3945974109459385552d-9, & -5.4853994841172088787d-8, 3.1601858228022739196d-6, & -3.3708782339998302320d-4, 1.6180215937704286491d-1 / data (d(i), i = 24, 31) / & -2.5640663123518180635d-14, 1.3288079339404032671d-12, & -4.3368537955908371563d-11, 1.3848103653102203186d-9, & -5.4824335664256344123d-8, 3.1601315173126153586d-6, & -3.3708776779035695640d-4, 1.6180215935248373474d-1 / data (d(i), i = 32, 39) / & -1.8678321325292127767d-14, 1.1354310934105733311d-12, & -4.1057197297998608931d-11, 1.3693990961296350970d-9, & -5.4762428935047089835d-8, 3.1599817092775027963d-6, & -3.3708756559715893599d-4, 1.6180215923508144240d-1 / if (x .lt. 0.86d0) then t = x * x y = ((((((a(0) * t + a(1)) * t + & a(2)) * t + a(3)) * t + a(4)) * t + & a(5)) * t + a(6)) * t + a(7) y = ((((((a(8) * t + a(9)) * t + & a(10)) * t + a(11)) * t + a(12)) * t + & a(13)) * t + a(14)) * t + a(15) - y * log(x) else if (x .lt. 4.15d0) then t = x - 5 / x k = int(t + 5) t = (k - 4) - t k = k * 14 y = ((((((((((((b(k) * t + b(k + 1)) * t + & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & b(k + 11)) * t + b(k + 12)) * t + b(k + 13) else if (x .lt. 12.5d0) then k = int(x) t = (k + 1) - x k = 15 * (k - 4) y = (((((((((((((c(k) * t + c(k + 1)) * t + & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & c(k + 11)) * t + c(k + 12)) * t + c(k + 13)) * t + & c(k + 14) else t = 60 / x k = 8 * (int(t)) y = (((((((d(k) * t + d(k + 1)) * t + & d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + & d(k + 5)) * t + d(k + 6)) * t + d(k + 7)) * & sqrt(t) * exp(-x) end if dbesk0 = y end Function dbesk0 ! ! Bessel K_1(x) function in double precision ! function dbesk1(x) implicit double precision (a - h, o - z) dimension a(0 : 15), b(0 : 119), c(0 : 119), d(0 : 39) data (a(i), i = 0, 15) / & 1.5151605362537935201d-13, 3.3637909513536510350d-11, & 5.6514041131016827202d-9, 6.7816840255069534052d-7, & 5.4253472222259226487d-5, 2.6041666666666637057d-3, & 6.2500000000000000090d-2, 5.0000000000000000000d-1, & -8.9790303384748696588d-11, -1.4029047449249185771d-8, & -1.5592893752540998113d-6, -1.1253607018469017569d-4, & -4.6421827665011579173d-3, -8.5370719728648409609d-2, & -3.0796575782920629660d-1, 1.0000000000000000004d0 / data (b(i), i = 0, 14) / & -9.4055461896630579928d-12, 3.1307934665844664773d-11, & 4.2005295001519243251d-10, -4.1636196779679820012d-9, & 1.4483026181700966164d-8, 1.1661000205428816914d-8, & -3.5023996724943046209d-7, 1.4404279316339005012d-6, & 5.3581564157158242080d-7, -3.5249754038612334639d-5, & 1.7150324075631641453d-4, -4.1276362611239191024d-5, & -4.6943110979636602591d-3, 3.5085369853392357659d-2, & 2.0063574339907819159d-1 / data (b(i), i = 15, 29) / & 3.3998989888944034586d-11, 7.1558979072937373055d-11, & -2.9226856932927698732d-9, 1.4591620256525610213d-8, & -6.6141635609854161666d-9, -1.9991101838984472332d-7, & 5.9185836628873530572d-7, 1.9562880347358085687d-6, & -1.5814366450418102764d-5, 7.6791682910944612028d-6, & 2.8354678948323983936d-4, -1.0217932669418690641d-3, & -3.2205661865210048433d-3, 4.3497494842354644077d-2, & 1.6110284302315089935d-1 / data (b(i), i = 30, 44) / & -2.0933987679665541827d-10, 7.9503322090520447316d-10,& 3.8000150948242575774d-9, -2.3076136195585571309d-8, & -2.3902880302550799653d-8, 3.1914500937804377478d-7, & 3.2639909831082417694d-7, -5.3166994792995439449d-6, & -3.1109524694269240094d-6, 9.2575906966353273247d-5, & 7.5063709094147644746d-7, -1.7416491592625765379d-3, & 1.2138560335171676007d-3, 4.5879687144659643175d-2, & 1.1566544716132846709d-1 / data (b(i), i = 45, 59) / & 3.1582384905164908749d-10, -1.9959561818098999516d-9, & 8.6959328920030927557d-10, 1.1642778282445577109d-8, & 4.3552264337818440471d-8, -1.5057982160481803238d-7, & -1.0101729117980989857d-6, 7.7002002510177612013d-7, & 1.9580574235590194233d-5, 1.9358461980242834361d-5, & -3.3932339942485532728d-4, -9.3416673584325090073d-4,& 5.5800080455912847227d-3, 3.8668683062477179235d-2, & 7.2651643500517000658d-2 / data (b(i), i = 60, 74) / & -1.1554749629758510059d-10, 8.2270678758893273006d-10, & -5.0211156951551538591d-10, -1.4929179050554858361d-9, & -2.7107940791526366702d-8, -4.2204764086705349384d-8, & 3.7253098167927628867d-7, 2.4374697215363361156d-6, & 1.4141942006909768370d-6, -4.8766389019473918231d-5, & -2.1681387247526720978d-4, 2.9325729929653405236d-4, & 6.4087534504827239815d-3, 2.6054628289709454356d-2, & 4.0156431128194184336d-2 / data (b(i), i = 75, 89) / & 2.5506555170746221691d-11, -1.3521164018407978152d-10, & -8.3281235274106699399d-11, -9.7764849575562351891d-10, & 3.4661828409940354542d-9, 3.9760633711791357544d-8, & 1.5902906645504529930d-7, -1.4919441249454941275d-7, & -5.3779684992094351263d-6, -2.7513862296246223142d-5, & -9.7880089725297162007d-6, 7.0787668964515789714d-4, & 4.6968199862345387583d-3, 1.4745740181663320127d-2, & 2.0048622219583455723d-2 / data (b(i), i = 90, 104) / & -3.4824483072529265585d-12, 1.5157161810563380451d-12, & 8.5303859696700686144d-12, 3.3455414203743741076d-10, & 2.0226016353844285376d-9, 5.3128154003266334990d-9, & -3.0799322316418042137d-8, -4.4455408272954712128d-7, & -2.4293274626893384034d-6, -3.2129079340119038354d-6, & 5.9225403683075388850d-5, 5.6822962576781683532d-4, & 2.7152446516406682732d-3, 7.4075873691848838485d-3, & 9.3044450815739269849d-3 / data (b(i), i = 105, 119) / & -2.7683216166276377232d-13, 3.1986676777610155465d-12, & 9.4142986954031445666d-12, 6.7934609179456399334d-11, & 3.4851529411470029330d-11, -2.5785248508896551557d-9, & -2.8310220027112571258d-8, -1.6384131113072271115d-7, & -3.2521663350596379097d-7, 4.0381388757622307160d-6, & 5.1917606978077281001d-5, 3.3420027947470126154d-4, & 1.3699550623118247094d-3, 3.4405619148342271096d-3, & 4.1042919106665762794d-3 / data (c(i), i = 0, 14) / & 4.5281968025889407937d-12, 1.0806749918195271176d-11, & 9.6200972728717669027d-11, 5.7214227063625263650d-10, & 3.6077804282954825099d-9, 2.2465236858536681852d-8, & 1.3676961264308735230d-7, 7.9561767489531997361d-7, & 4.3014380065615550573d-6, 2.0921713905550285590d-5, & 8.8079183950590176926d-5, 3.0549414408830252064d-4, & 8.1295715613927890473d-4, 1.4679809476357079195d-3, & 1.3439197177355090057d-3 / data (c(i), i = 15, 29) / & 7.6019964430402432637d-13, -2.2616198599158271190d-13, & 1.7904450823779000744d-11, 9.1467054855312232717d-11, & 7.1378582044879519122d-10, 4.9925255415445769102d-9, & 3.3767315471315546644d-8, 2.1350774539167751457d-7, & 1.2314353082655232903d-6, 6.2918685053670619181d-6, & 2.7493229298777000013d-5, 9.8085825401369821771d-5, & 2.6670282677770444935d-4, 4.8967895428135985381d-4, & 4.5418248688489697144d-4 / data (c(i), i = 30, 44) / & 9.4180115230375147213d-14, 7.5943117003734061145d-14, & 3.0335730243874287654d-12, 2.0202796115462268051d-11, & 1.6839020189186971198d-10, 1.2907875663127201526d-9, & 9.3547676125865798920d-9, 6.2471974110281880722d-8, & 3.7585985422997380441d-7, 1.9838348288114906484d-6, & 8.8884862203671982034d-6, 3.2333259238682810218d-5, & 8.9266668913380400243d-5, 1.6589185669844051903d-4, & 1.5536921180500113394d-4 / data (c(i), i = 45, 59) / & 1.5425475332301107271d-14, 2.8674534590132490434d-14, & 6.5078462279160216936d-13, 5.0939757793961391211d-12, & 4.4979837460748975520d-11, 3.6662925847520171711d-10, & 2.7848878755089582413d-9, 1.9298120059339477820d-8, & 1.1950323861976892013d-7, 6.4513432758147478287d-7, & 2.9422095033982461936d-6, 1.0854433321174584937d-5, & 3.0307433185818899481d-5, 5.6840981443065017850d-5, & 5.3637016379451945253d-5 / data (c(i), i = 60, 74) / & 3.1077953698439839352d-15, 8.6899496170729520378d-15, & 1.6258562067326054104d-13, 1.4104842571366761537d-12, & 1.3019455544084110747d-11, 1.1070466372863950239d-10, & 8.6890603844230597917d-10, 6.1793722175049967488d-9, & 3.9058865943755615801d-8, 2.1432806981070368523d-7, & 9.9034657762983230155d-7, 3.6925185861895664251d-6, & 1.0399877577259449786d-5, 1.9644939661550210015d-5, & 1.8648773453825584597d-5 / data (c(i), i = 75, 89) / & 7.2831555285336286457d-16, 2.6077534095895783532d-15, & 4.4881202059263153495d-14, 4.1674329383944385626d-13, & 3.9760100480223728037d-12, 3.4835976355351183010d-11, & 2.7993254212770249700d-10, 2.0286513276830758107d-9, & 1.3018343087118439152d-8, 7.2315927974997999365d-8, & 3.3750708681924201599d-7, 1.2688020879407355571d-6, & 3.5980954090811587848d-6, 6.8358260635246667316d-6, & 6.5208606745808860557d-6 / data (c(i), i = 90, 104) / & 1.9026412343503745875d-16, 8.0073765508732553766d-16, & 1.3245754278055523992d-14, 1.2885201653055058502d-13, & 1.2600129301230402587d-12, 1.1283306843147549277d-11, & 9.2261481309646814329d-11, 6.7812033168299846818d-10, & 4.4020645304595102132d-9, 2.4685719238301517679d-8, & 1.1611886719473112509d-7, 4.3940380936523135466d-7, & 1.2529878285546791905d-6, 2.3917218527087570384d-6, & 2.2907574647671878160d-6 / data (c(i), i = 105, 119) / & 5.3709522135744366512d-17, 2.5239221050372845433d-16, & 4.0933147145899083360d-15, 4.1152784247617592367d-14, & 4.0998840572769381012d-13, 3.7319354625807158852d-12, & 3.0921671702920868014d-11, 2.2975898538634445343d-10, & 1.5049754445782364328d-9, 8.5030864719789148982d-9, & 4.0250559391118423810d-8, 1.5312755642491878591d-7, & 4.3865020375297892208d-7, 8.4059737392822153101d-7, & 8.0785884122023473319d-7 / data (d(i), i = 0, 7) / & 9.2371554649979581914d-14, -2.3111336195788410887d-12, & 5.7728710326649832559d-11, -1.8002298130091372598d-9, & 7.6810375010517145638d-8, -5.2669973752193823306d-6, & 1.0112634961227401357d-3, 1.6180215937964160466d-1 / data (d(i), i = 8, 15) / & 6.1381146507252683381d-14, -2.1034499679806301862d-12, & 5.7090233460448415278d-11, -1.7990724350642330817d-9, & 7.6809056078388019946d-8, -5.2669964425290062357d-6, & 1.0112634957478283390d-3, 1.6180215937970716383d-1 / data (d(i), i = 16, 23) / & 4.2458150578401296419d-14, -1.8435733128339016981d-12, & 5.5534955081564656595d-11, -1.7938162188526358466d-9, & 7.6798230945934117807d-8, -5.2669828728791921259d-6, & 1.0112634861753356559d-3, 1.6180215938263409582d-1 / data (d(i), i = 24, 31) / & 3.0314798962267007518d-14, -1.5915009905364214455d-12, & 5.3275907427402047438d-11, -1.7824862013841369751d-9, & 7.6763890356447075810d-8, -5.2669199860465945909d-6, & 1.0112634217687349189d-3, 1.6180215941108227283d-1 / data (d(i), i = 32, 39) / & 2.2211515002229271212d-14, -1.3664088221521734796d-12, & 5.0585177270502341602d-11, -1.7645432205894533462d-9, & 7.6691805594577373698d-8, -5.2667455286976269634d-6, & 1.0112631862810974580d-3, 1.6180215954783127877d-1 / if (x .lt. 0.8d0) then t = x * x y = (((((((a(0) * t + a(1)) * t + & a(2)) * t + a(3)) * t + a(4)) * t + & a(5)) * t + a(6)) * t + a(7)) * x y = (((((((a(8) * t + a(9)) * t + & a(10)) * t + a(11)) * t + a(12)) * t + & a(13)) * t + a(14)) * t + a(15)) / x + & y * log(x) else if (x .lt. 5.5d0) then v = 3 / x t = x - v k = int(t + 3) t = (k - 2) - t k = k * 15 y = ((((((((((((((b(k) * t + b(k + 1)) * t + & b(k + 2)) * t + b(k + 3)) * t + b(k + 4)) * t + & b(k + 5)) * t + b(k + 6)) * t + b(k + 7)) * t + & b(k + 8)) * t + b(k + 9)) * t + b(k + 10)) * t + & b(k + 11)) * t + b(k + 12)) * t + b(k + 13)) * t + & b(k + 14)) * v else if (x .lt. 12.5d0) then k = int(x) t = (k + 1) - x k = 15 * (k - 5) y = (((((((((((((c(k) * t + c(k + 1)) * t + & c(k + 2)) * t + c(k + 3)) * t + c(k + 4)) * t + & c(k + 5)) * t + c(k + 6)) * t + c(k + 7)) * t + & c(k + 8)) * t + c(k + 9)) * t + c(k + 10)) * t + & c(k + 11)) * t + c(k + 12)) * t + c(k + 13)) * t + & c(k + 14) else t = 60 / x k = 8 * (int(t)) y = (((((((d(k) * t + d(k + 1)) * t + & d(k + 2)) * t + d(k + 3)) * t + d(k + 4)) * t + & d(k + 5)) * t + d(k + 6)) * t + d(k + 7)) * & sqrt(t) * exp(-x) end if dbesk1 = y end Function dbesk1 ! End Module BESSEL !********************************************************************** !********************************************************************** Module ZERO ! Implicit None ! PRIVATE PUBLIC :: BRENT Integer, Parameter :: DP = KIND(1.d0) ! CONTAINS ! Real (DP) Function BRENT ( a, b, machep, t, F ) ! Trouvé sur http://people.sc.fsu.edu/~burkardt/f_src/brent/brent.f90 ! Real(DP), Intent(In) :: a, b, machep, t ! Interface Double Precision Function F( X ) Result ( y ) Double Precision, Intent(In) :: X End Function F End Interface ! !*****************************************************************************80 ! !! ZERO seeks the root of a function F(X) in an interval [A,B]. ! ! Discussion: ! ! The interval [A,B] must be a change of sign interval for F. ! That is, F(A) and F(B) must be of opposite signs. Then ! assuming that F is continuous implies the existence of at least ! one value C between A and B for which F(C) = 0. ! ! The location of the zero is determined to within an accuracy ! of 6 * MACHEPS * abs ( C ) + 2 * T. ! ! Licensing: ! ! This code is distributed under the GNU LGPL license. ! ! Modified: ! ! 12 April 2008 ! ! Author: ! ! Original FORTRAN77 version by Richard Brent ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Richard Brent, ! Algorithms for Minimization Without Derivatives, ! Dover, 2002, ! ISBN: 0-486-41998-3, ! LC: QA402.5.B74. ! ! Parameters: ! ! Input, real ( kind = 8 ) A, B, the endpoints of the change of sign interval. ! ! Input, real ( kind = 8 ) MACHEP, an estimate for the relative machine ! precision. ! ! Input, real ( kind = 8 ) T, a positive error tolerance. ! ! Input, external real ( kind = 8 ) F, the name of a user-supplied ! function, of the form "FUNCTION F ( X )", which evaluates the ! function whose zero is being sought. ! ! Output, real ( kind = 8 ) ZERO, the estimated value of a zero of ! the function F. ! ! real ( kind = 8 ) c real ( kind = 8 ) d real ( kind = 8 ) e !real ( kind = 8 ) f real ( kind = 8 ) fa real ( kind = 8 ) fb real ( kind = 8 ) fc real ( kind = 8 ) m real ( kind = 8 ) p real ( kind = 8 ) q real ( kind = 8 ) r real ( kind = 8 ) s real ( kind = 8 ) sa real ( kind = 8 ) sb real ( kind = 8 ) tol ! ! Make local copies of A and B. ! sa = a sb = b fa = F ( sa ) fb = F ( sb ) c = sa fc = fa e = sb - sa d = e do if ( abs ( fc ) < abs ( fb ) ) then sa = sb sb = c c = sa fa = fb fb = fc fc = fa end if tol = 2.0D+00 * machep * abs ( sb ) + t m = 0.5D+00 * ( c - sb ) if ( abs ( m ) <= tol .or. fb == 0.0D+00 ) then exit end if if ( abs ( e ) < tol .or. abs ( fa ) <= abs ( fb ) ) then e = m d = e else s = fb / fa if ( sa == c ) then p = 2.0D+00 * m * s q = 1.0D+00 - s else q = fa / fc r = fb / fc p = s * ( 2.0D+00 * m * a * ( q - r ) - ( sb - sa ) * ( r - 1.0D+00 ) ) q = ( q - 1.0D+00 ) * ( r - 1.0D+00 ) * ( s - 1.0D+00 ) end if if ( 0.0D+00 < p ) then q = - q else p = - p end if s = e e = d if ( 2.0D+00 * p < 3.0D+00 * m * q - abs ( tol * q ) .and. & p < abs ( 0.5D+00 * s * q ) ) then d = p / q else e = m d = e end if end if sa = sb fa = fb if ( tol < abs ( d ) ) then sb = sb + d else if ( 0.0D+00 < m ) then sb = sb + tol else sb = sb - tol end if fb = F ( sb ) if ( ( 0.0D+00 < fb .and. 0.0D+00 < fc ) .or. & ( fb <= 0.0D+00 .and. fc <= 0.0D+00 ) ) then c = sa fc = fa e = sb - sa d = e end if end do BRENT = sb return end function BRENT ! End Module ZERO !********************************************************************** !********************************************************************** Module INTEGRALE ! Implicit None ! PRIVATE PUBLIC :: SIMPSON, GAUSS_LAGUERRE_8PTS, GAUSS_LAGUERRE_12PTS ! Integer, Parameter :: SP = KIND(1.0) , DP = KIND(1.d0) , P = DP ! CONTAINS !------------------------------------------------------------------- Real(P) Function SIMPSON ( F , a , b , n ) Result ( y ) ! Interface Function F ( x ) Result ( y ) Integer, Parameter :: SP = KIND(1.0) , DP = KIND(1.d0) , P = DP ! pas encore en F2003 Real(P), Intent(In) :: x Real(P) :: y End Function F End Interface ! Real(P), Intent(In) :: a, b Integer, Intent(In) :: n ! Real(P) :: h, x_i, y_i, y_a, y_b Integer :: i ! h = 0.5_P * (b-a) / REAL(n,P) y_a = F(a) y_b = F(b) ! If ( y_a < y_b ) Then y = y_a Do i = 1 , 2*n-1 x_i = a + h * REAL(i,P) y_i = F(x_i) y = y + 2.0_P * REAL( 1 + MODULO(i,2) , P ) * y_i End Do y = h * ( y + y_b ) / 3.0_P Else y = y_b Do i = 2*n-1 , 1 , -1 x_i = b - h * REAL(i,P) y_i = F(x_i) y = y + 2.0_P * REAL( 1 + MODULO(i,2) , P ) * y_i End Do y = h * ( y + y_a ) / 3.0_P End If ! End Function SIMPSON !------------------------------------------------------------------- Real(P) Function GAUSS_LAGUERRE_8PTS( F ) Result(y) ! Approximate the integral of F(x) exp(-x) from 0 to infinity using the ! 8 point Gauss-Laguerre integral approximation formula. ! Interface Function F ( x ) Result ( y ) Integer, Parameter :: SP = KIND(1.0) , DP = KIND(1.d0) , P = DP ! pas encore en F2003 Real(P), Intent(In) :: x Real(P) :: y End Function F End Interface ! Integer :: i ! ! Coefficients pris sur http://mymathlib.webtrellis.net/c_source/quadrature/gauss_laguerre/gauss_laguerre_8pts.c Real(P), dimension(8), parameter :: x_i = (/ & 1.70279632305100999786D-01, 9.03701776799379912170D-01, & 2.25108662986613068929D+00, 4.26670017028765879378D+00, & 7.04590540239346569719D+00, 1.07585160101809952241D+01, & 1.57406786412780045781D+01, 2.28631317368892641052D+01 /) Real(P), dimension(8), parameter :: w_i = (/ & 3.69188589341637529929D-01, 4.18786780814342956078D-01, & 1.75794986637171805706D-01, 3.33434922612156515224D-02, & 2.79453623522567252491D-03, 9.07650877335821310457D-05, & 8.48574671627253154502D-07, 1.04800117487151038157D-09 /) ! y = 0 Do i = 1, 8 y = y + w_i(i) * F( x_i(i) ) End Do ! End Function GAUSS_LAGUERRE_8PTS !------------------------------------------------------------------- Real(P) Function GAUSS_LAGUERRE_12PTS( F ) Result(y) ! Approximate the integral of F(x) exp(-x) from 0 to infinity using the ! 12 point Gauss-Laguerre integral approximation formula. ! Interface Function F ( x ) Result ( y ) Integer, Parameter :: SP = KIND(1.0) , DP = KIND(1.d0) , P = DP ! pas encore en F2003 Real(P), Intent(In) :: x Real(P) :: y End Function F End Interface ! Integer :: I ! ! Coefficients pris sur http://mymathlib.webtrellis.net/c_source/quadrature/gauss_laguerre/gauss_laguerre_12pts.c Real(P), dimension(12), parameter :: x_i = (/ & 1.15722117358020675267D-01, 6.11757484515130665371D-01, & 1.51261026977641878681D+00, 2.83375133774350722858D+00, & 4.59922763941834848462D+00, 6.84452545311517734786D+00, & 9.62131684245686704358D+00, 1.30060549933063477205D+01, & 1.71168551874622557277D+01, 2.21510903793970056700D+01, & 2.84879672509840003117D+01, 3.70991210444669203376D+01 /) Real(P), dimension(12), parameter :: w_i = (/ & 2.64731371055443190344D-01, 3.77759275873137982026D-01, & 2.44082011319877564251D-01, 9.04492222116809307247D-02, & 2.01023811546340965218D-02, 2.66397354186531588103D-03, & 2.03231592662999392117D-04, 8.36505585681979874511D-06, & 1.66849387654091026123D-07, 1.34239103051500414548D-09, & 3.06160163503502078133D-12, 8.14807746742624168225D-16 /) ! y = 0 Do I = 1, 12 y = y + w_i(I) * F( x_i(I) ) End Do ! End Function GAUSS_LAGUERRE_12PTS !------------------------------------------------------------------- End Module INTEGRALE