1 | SUBROUTINE PIACH4(NU,NT,XROT,XTRA,T) |
---|
2 | IMPLICIT REAL (A-H,O-Z) |
---|
3 | REAL NU,I10,NSQR,K |
---|
4 | DIMENSION T(100) |
---|
5 | DIMENSION FOX(60),XJ(23),G(23),E(23),WW(3,23),RHO(100,23), |
---|
6 | $I10(100),XROT(100),XTRA(100) |
---|
7 | |
---|
8 | save ww,a2,pi,c,hck,tau1,tau2,cg,rho,c1,k,asqr,quads1,i10,xj |
---|
9 | |
---|
10 | DATA FOX /1.00,1.00,4.12,1.00,0.52,0.81,1.00,0.70,0.75,0.80, |
---|
11 | $0.69,1.24,0.44,1.29,1.05,0.56,0.79,0.79,0.89,1.16,1.27,0.86, |
---|
12 | $0.95,1.00,0.72,1.16,0.92,1.00,0.92,1.18,0.90,1.18,1.09,0.85, |
---|
13 | $0.96,0.93,0.82,1.18,1.17,0.95,0.96,1.00,0.85,1.06,0.90,15*1.00/ |
---|
14 | DATA JMAX1/20/,JMAX2/23/,PI/3.141592654/,H/1.05459E-27/, |
---|
15 | $C/2.997925E10/,K /1.38054E-16/,HCK/1.43892/,XNLOS/2.687E19/, 276. |
---|
16 | $ASQR/6.7081/,EPSK/148.6/,B/5.25/,ANORM/4.797E17/ |
---|
17 | DO 8 IT=1,NT |
---|
18 | X=4.*EPSK/T(IT) |
---|
19 | I10(IT)=FI10(X) |
---|
20 | 8 CONTINUE |
---|
21 | DO 10 J=1,JMAX2 |
---|
22 | XJ(J)=FLOAT(J-1) |
---|
23 | E(J)=B*XJ(J)*(XJ(J)+1.) |
---|
24 | 10 CONTINUE |
---|
25 | G(1)=5. |
---|
26 | G(2)=3. |
---|
27 | G(3)=5. |
---|
28 | G(4)=11. |
---|
29 | G(5)=13. |
---|
30 | G(6)=11. |
---|
31 | !!! correction au nez |
---|
32 | DO 15 J=7,JMAX2 |
---|
33 | 15 G(J)=G(J-6)+16. |
---|
34 | DO 20 J=1,JMAX1 |
---|
35 | DO 20 I=1,3 |
---|
36 | WW(I,J)=2.*PI*C*(E(J+I)-E(J)) |
---|
37 | 20 CONTINUE |
---|
38 | NSQR=(XNLOS*1.0E-24)**2 |
---|
39 | c C1=64./63.*PI*PI*NSQR/H/C/ANORM |
---|
40 | C1=64./63.*PI*PI*NSQR/H/C |
---|
41 | QUADS1=11.10 |
---|
42 | TAU1=9.00E-14 |
---|
43 | TAU2=13.60E-14 |
---|
44 | DO 90 IT=1,NT |
---|
45 | SUM=0. |
---|
46 | DO 50 J=1,JMAX2 |
---|
47 | ARG=HCK*E(J)/T(IT) |
---|
48 | RHO(IT,J)=EXP(-ARG) |
---|
49 | SUM=SUM+(2.*XJ(J)+1.)*G(J)*RHO(IT,J) |
---|
50 | 50 CONTINUE |
---|
51 | DO 55 J=1,JMAX2 |
---|
52 | RHO(IT,J)=RHO(IT,J)/SUM |
---|
53 | 55 CONTINUE |
---|
54 | 90 CONTINUE |
---|
55 | RETURN |
---|
56 | ENTRY OPACH4(NU,NT,XROT,XTRA,T) |
---|
57 | W=2.*PI*C*NU |
---|
58 | DO 200 IT=1,NT |
---|
59 | XROT(IT)=0. |
---|
60 | XTRA(IT)=0. |
---|
61 | CZ=HCK*NU/T(IT) |
---|
62 | CW=W*(1.-EXP(-CZ)) |
---|
63 | DO 125 J=1,JMAX1 |
---|
64 | XTRA(IT)=XTRA(IT)+(2.*XJ(J)+1.)*(2.*XJ(J)+1.)* |
---|
65 | $ RHO(IT,J)*GAMMB(W,T(IT),TAU1,TAU2) |
---|
66 | !!! avant cette boucle contenait K au lieu de IK en contradiction |
---|
67 | !!! avec K Boltzman |
---|
68 | DO 25 I=1,3 |
---|
69 | IK=3*(J-1)+I |
---|
70 | XROT(IT)=XROT(IT)+FOX(IK)*(2.*XJ(J)+1.)*(2.*XJ(J+I)+1.)* |
---|
71 | $(RHO(IT,J )*GAMMB(W-WW(I,J),T(IT),TAU1,TAU2) |
---|
72 | $+RHO(IT,J+I)*GAMMB(W+WW(I,J),T(IT),TAU1,TAU2)) |
---|
73 | 25 CONTINUE |
---|
74 | 125 CONTINUE |
---|
75 | c COEF=CW*C1*IK*ASQR*QUADS1*I10(IT) |
---|
76 | COEF=CW*C1*K*ASQR*QUADS1*I10(IT) |
---|
77 | XROT(IT)=XROT(IT)*COEF |
---|
78 | XTRA(IT)=XTRA(IT)*COEF |
---|
79 | 200 CONTINUE |
---|
80 | RETURN |
---|
81 | END |
---|