source: trunk/LMDZ.TITAN.old/libf/phytitan/gamfcn.F @ 3094

Last change on this file since 3094 was 3, checked in by slebonnois, 14 years ago

Creation de repertoires:

  • chantiers : pour communiquer sur nos projets de modifs
  • documentation : pour stocker les docs

Ajout de:

  • libf/phytitan : physique de Titan
  • libf/chimtitan: chimie de Titan
  • libf/phyvenus : physique de Venus
File size: 2.6 KB
Line 
1      FUNCTION GAMFCN(W,T,TAU1,TAU2)
2      IMPLICIT REAL (A-H,O-Z)
3C
4C     COMPUTES THE LINE SHAPE FOR PRESSURE-INDUCED H2-H2 AND H2-HE
5C     TRANSITIONS, FROM THE SEMI-EMIRICAL FORMULAE OF BIRNBAUM ET AL.
6C     SEE EG., GEORGE BIRNBAUM AND E. RICHARD COHEN, CANADIAN JOURNAL
7C     OF PHYSICS, VOL. 54, 593 (1976).
8C     NOTE THAT "GAMFCN" IS A POOR NAME FOR THIS ROUTINE; THE GAMMA
9C     OF BIRNBAUM AND COHEN IS NOT THE USUAL "GAMMA FUNCTION".
10C
11      save hk,pi,tau10,tau20
12      DATA HK/7.638315E-12/,PI/3.141593/,TAU10/0.0/,TAU20/0.0/
13      logical first
14      data first/.true./
15      save first
16
17C     NOTE: HK = 1.05450E-27 / 1.38054E-16
18C
19C***********************************************************************
20      save tau12,z2,hbh
21
22      if (first)
23     s  print*,'WARNING!!! ON rajoute des valeurs a 0.',
24     s   'Est-ce bien raisonable... dans GAMFCN'
25      first=.false.
26      IF (TAU1 .NE. TAU10) GO TO 10
27      IF (TAU2 .EQ. TAU20) GO TO 20
28 10   TAU12 = TAU1 * TAU1
29      TAU22 = TAU2 * TAU2
30      HBH = 0.5 * HK / T
31      Z2 = SQRT(TAU22 + HBH**2) / TAU1
32      TAU10 = TAU1
33      TAU20 = TAU2
34 20   WSQR = W * W
35      Z = SQRT(1.0+WSQR*TAU12) * Z2
36      IF (Z .LE. 1.0) GO TO 50
37C     COMPUTE K1 BESSEL FUNCTION USING POLYNOMIAL APPROXIMATION
38      A = 1.0 / Z
39      BK1 = 1.253314 + .4699927*A
40      B = A * A
41      BK1 = BK1 - .1468583*B
42      B = B * A
43      BK1 = BK1 + .1280427*B
44      B = B * A
45      BK1 = BK1 - .1736432*B
46      B = B * A
47      BK1 = BK1 + .2847618*B
48      B = B * A
49      BK1 = BK1 - .4594342*B
50      B = B * A
51      BK1 = BK1 + .6283381*B
52      B = B * A
53      BK1 = BK1 - .6632295*B
54      B = B * A
55      BK1 = BK1 + .5050239*B
56      B = B * A
57      BK1 = BK1 - .2581304*B
58      B = B * A
59      BK1 = BK1 + .7880001E-01*B
60      B = B * A
61      BK1 = BK1 - .1082418E-01*B
62      BK1 = EXP(-Z) * BK1 * SQRT(A)
63      GO TO 100
64C     COMPUTE K1 BESSEL FUNCTION USING SERIES EXPANSION
65 50   A = 0.5 * Z
66      B = .5772157 + LOG(A)
67      C = A * A
68      BK1 = 1.0/Z + A*(B-0.5)
69      A = A * C
70      BK1 = BK1 + A*.2500000E+00*(0.5+(B-1.500000)*2.0)
71      A = A * C
72      BK1 = BK1 + A*.2777777E-01*(0.5+(B-1.833333)*3.0)
73      A = A * C
74      BK1 = BK1 + A*.1736110E-02*(0.5+(B-2.083333)*4.0)
75      A = A * C
76      BK1 = BK1 + A*.6944439E-04*(0.5+(B-2.283333)*5.0)
77      A = A * C
78      BK1 = BK1 + A*.1929009E-05*(0.5+(B-2.449999)*6.0)
79      A = A * C
80      BK1 = BK1 + A*.3936752E-07*(0.5+(B-2.592855)*7.0)
81      A = A * C
82      BK1 = BK1 + A*.6151173E-09*(0.5+(B-2.717855)*8.0)
83 100  CONTINUE
84      GAMFCN = TAU1/PI * EXP(TAU2/TAU1+HBH*W) * Z*BK1 / (1.0+WSQR*TAU12)
85      RETURN
86      END
Note: See TracBrowser for help on using the repository browser.