source: trunk/fsource/pypowder.for @ 5369

Last change on this file since 5369 was 5369, checked in by toby, 7 months ago

Incorporate James Hester's updates for high angle FCJ
File size: 8.9 KB
Line 
1      Subroutine PYPSVFCJ(NPTS,DTT,TTHETA,CDSIG,CDGAM,SPH,PRFUNC)
2C DTT in degrees
3C TTHETA in degrees
4C     SPH is S/L + H/L
5C     CDSIG Gaussian variance,centidegrees squared
6C     CDGAM Lorenzian FWHM, centidegrees
7C RETURNS FUNCTION ONLY
8Cf2py intent(in) NPTS
9Cf2py intent(in) DTT
10Cf2py depend(NPTS) DTT
11Cf2py intent(in) TTHETA
12Cf2py intent(in) CDSIG
13Cf2py intent(in) CDGAM
14Cf2py intent(in) SPH
15Cf2py intent(out) PRFUNC
16Cf2py depend(NPTS) PRFUNC
17      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
18      REAL*4 TTHETA,SIG,CDSIG,CDGAM,SPH,DPRDT,DPRDG,DPRDD,DPRDLZ,DPRDS
19      REAL*4 GAM
20      INTEGER*4 NPTS,I
21C     CDSIG is in centidegrees squared, we must change to normal degrees
22      SIG = CDSIG/10000.0
23      GAM = CDGAM/100.0
24      DO I=0,NPTS-1
25         Call Get_Prof_Val(SIG,GAM,SPH,0.0,DTT(I)+TTHETA,TTHETA,DPRDT,
26     1    DPRDG,DPRDD,DPRDS,DPRDLZ,PRFUNC(I))
27          PRFUNC(I)=PRFUNC(I)/100.   !Calling code expects peak normalised in centidegrees
28      END DO
29      RETURN
30      END
31
32      SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,CDSIG,CDGAM,SPH,PRFUNC,
33     1   DPRDT,SIGPART,GAMPART,SLPART)
34C DTT in degrees
35C TTHETA in degrees
36C SPH is S/L + H/L
37C RETURNS FUNCTION & DERIVATIVES
38Cf2py intent(in) NPTS
39Cf2py intent(in) DTT
40Cf2py depend(NPTS) DTT
41Cf2py intent(in) TTHETA
42Cf2py intent(in) CDSIG
43Cf2py intent(in) CDGAM
44Cf2py intent(in) SPH
45Cf2py intent(out) PRFUNC
46Cf2py depend(NPTS) PRFUNC
47Cf2py intent(out) DPRDT
48Cf2py depend(NPTS) DPRDT
49Cf2py intent(out) SIGPART
50Cf2py depend(NPTS) SIGPART
51Cf2py intent(out) GAMPART
52Cf2py depend(NPTS) GAMPART
53Cf2py intent(out) SLPART
54Cf2py depend(NPTS) SLPART
55      INTEGER*4 NPTS,I
56      REAL*4 TTHETA,CDSIG,SIG,CDGAM,SPH,LPART
57      REAL*4 GAM
58      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1)
59      REAL*4 GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
60      SIG = CDSIG/10000.
61      GAM = CDGAM/100.
62      DO I=0,NPTS-1
63         Call Get_Prof_Val(SIG,GAM,SPH,0.0,DTT(I)+TTHETA,TTHETA,
64     1      DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I),LPART, PRFUNC(I))
65         ! Calling code expects all values to be for a peak normalised in centidegrees
66         SIGPART(I)=SIGPART(I)/1.0e6
67         GAMPART(I)=GAMPART(I)/1.0e4
68         SLPART(I)=SLPART(I)/100.
69         PRFUNC(I)=PRFUNC(I)/100.
70         DPRDT(I)=DPRDT(I)/100.
71      END DO
72      RETURN
73      END
74
75      SUBROUTINE PYPSVOIGT(NPTS,DTT,SIG,GAM,PRFUNC)
76C DTT in degrees
77C RETURNS FUNCTION ONLY
78Cf2py intent(in) NPTS
79Cf2py intent(in) DTT
80Cf2py depend(NPTS) DTT
81Cf2py intent(in) SIG
82Cf2py intent(in) GAM
83Cf2py intent(out) PRFUNC
84Cf2py depend(NPTS) PRFUNC
85
86      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
87      REAL*4 SIG,GAM
88      INTEGER*4 NPTS,I
89      DO I=0,NPTS-1
90        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
91     1     PRFUNC(I),DPRDT,SIGPART,GAMPART)
92      END DO
93      RETURN
94      END
95
96      SUBROUTINE PYDPSVOIGT(NPTS,DTT,SIG,GAM,PRFUNC,
97     1   DPRDT,SIGPART,GAMPART)
98C DTT in degrees
99C RETURNS FUNCTION & DERIVATIVES
100Cf2py intent(in) NPTS
101Cf2py intent(in) DTT
102Cf2py depend(NPTS) DTT
103Cf2py intent(in) SIG
104Cf2py intent(in) GAM
105Cf2py intent(out) PRFUNC
106Cf2py depend(NPTS) PRFUNC
107Cf2py intent(out) DPRDT
108Cf2py depend(NPTS) DPRDT
109Cf2py intent(out) SIGPART
110Cf2py depend(NPTS) SIGPART
111Cf2py intent(out) GAMPART
112Cf2py depend(NPTS) GAMPART
113
114      INTEGER*4 NPTS
115      REAL*4 SIG,GAM
116      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
117     1   GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
118      DO I=0,NPTS-1
119        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
120     1     PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I))
121        DPRDT(I) = DPRDT(I)*100.
122      END DO
123      RETURN
124      END
125
126      SUBROUTINE PYPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
127C DTT in degrees
128C TTHETA in degrees
129C SPH is S/L + H/L
130C RETURNS FUNCTION ONLY
131Cf2py intent(in) NPTS
132Cf2py intent(in) DTT
133Cf2py depend(NPTS) DTT
134Cf2py intent(in) TTHETA
135Cf2py intent(in) SIG
136Cf2py intent(in) GAM
137Cf2py intent(in) SPH
138Cf2py intent(out) PRFUNC
139Cf2py depend(NPTS) PRFUNC
140
141      INTEGER*4 NPTS
142      REAL*4 TTHETA,SIG,GAM,SPH
143      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
144      DO I=0,NPTS-1
145        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
146     1     PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
147      END DO
148      RETURN
149      END
150
151      SUBROUTINE PYDPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SHL,PRFUNC,
152     1  DPRDT,SIGPART,GAMPART,SLPART)
153C DTT in degrees
154C TTHETA in degrees
155C SPH is S/L + H/L
156C RETURNS FUNCTION & DERIVATIVES
157Cf2py intent(in) NPTS
158Cf2py intent(in) DTT
159Cf2py depend(NPTS) DTT
160Cf2py intent(in) TTHETA
161Cf2py intent(in) SIG
162Cf2py intent(in) GAM
163Cf2py intent(in) SHL
164Cf2py intent(out) PRFUNC
165Cf2py depend(NPTS) PRFUNC
166Cf2py intent(out) DPRDT
167Cf2py depend(NPTS) DPRDT
168Cf2py intent(out) SIGPART
169Cf2py depend(NPTS) SIGPART
170Cf2py intent(out) GAMPART
171Cf2py depend(NPTS) GAMPART
172Cf2py intent(out) SLPART
173Cf2py depend(NPTS) SLPART
174
175      INTEGER*4 NPTS
176      REAL*4 TTHETA,SIG,GAM,SHL
177      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
178     1   GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
179      DO I=0,NPTS-1
180        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
181     1     PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
182          SLPART(I) = SPART
183        DPRDT(I) = DPRDT(I)*100.
184      END DO
185      RETURN
186      END
187
188      SUBROUTINE PYEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC)
189C DTT in microsec
190C RETURNS FUNCTION
191Cf2py intent(in) NPTS
192Cf2py intent(in) DTT
193Cf2py depend(NPTS) DTT
194Cf2py intent(in) ALP
195Cf2py intent(in) BET
196Cf2py intent(in) SIG
197Cf2py intent(in) GAM
198Cf2py intent(out) PRFUNC
199Cf2py depend(NPTS) PRFUNC
200
201      INTEGER*4 NPTS
202      REAL*4 ALP,BET,SIG,GAM,SHL
203      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1),DPRDT,ALPPART,
204     1   BETPART,SIGPART,GAMPART
205      DO I=0,NPTS-1
206        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT,
207     1     ALPPART,BETPART,SIGPART,GAMPART)
208      END DO
209      RETURN
210      END
211
212      SUBROUTINE PYDEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC,
213     1   DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
214C DTT in microsec
215C RETURNS FUNCTION & DERIVATIVES
216Cf2py intent(in) NPTS
217Cf2py intent(in) DTT
218Cf2py depend(NPTS) DTT
219Cf2py intent(in) ALP
220Cf2py intent(in) BET
221Cf2py intent(in) SIG
222Cf2py intent(in) GAM
223Cf2py intent(out) PRFUNC
224Cf2py depend(NPTS) PRFUNC
225Cf2py intent(out) DPRDT
226Cf2py depend(NPTS) DPRDT
227Cf2py intent(out) ALPPART
228Cf2py depend(NPTS) ALPPART
229Cf2py intent(out) BETPART
230Cf2py depend(NPTS) BETPART
231Cf2py intent(out) SIGPART
232Cf2py depend(NPTS) SIGPART
233Cf2py intent(out) GAMPART
234Cf2py depend(NPTS) GAMPART
235
236      INTEGER*4 NPTS
237      REAL*4 ALP,BET,SIG,GAM,SHL
238      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),ALPPART(0:NPTS-1),
239     1   BETPART(0:NPTS-1),SIGPART(0:NPTS-1), 
240     1   GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
241      DO I=0,NPTS-1
242        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT(I),
243     1     ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
244      END DO
245      RETURN
246      END
247     
248C python interface to gauleg
249      SUBROUTINE PYGAULEG(X1,X2,N,GL,WT)
250Cf2py intent(in) X1
251Cf2py intent(in) X2
252Cf2py intent(in) N
253Cf2py intent(out) GL
254Cf2py depend(N) GL
255Cf2py intent(out) WT
256Cf2py depends(N) WT
257
258      INTEGER*4 N
259      REAL*4 X1,X2,GL(0:N-1),WT(0:N-1)
260      CALL GAULEG(X1,X2,GL,WT,N)
261      RETURN
262      END
263
264
265C Fortran (fast) linear interpolation -- B.H. Toby 9/2011
266      SUBROUTINE PYFINTERP(NIN,XIN,YIN,NOUT,XOUT,YOUT)
267C XIN(1:NIN) and YIN(1:NIN) are arrays of (x,y) points to be interpolated
268C     Values must be sorted increasing in XIN
269C XOUT(1:NOUT) is an array of x values, must also be sorted increasing in x
270C    XOUT may contain values smaller than XIN(1) or larger than XIN(NIN)
271C RETURNS interpolated y values corresponding to XOUT. Values outside the
272C   range of XIN are set to zero.
273C Needs a way to signal an error if XIN or XOUT is not sorted -- for now stops
274Cf2py intent(in) NIN
275Cf2py intent(in)  XIN
276Cf2py depend(NIN) XIN
277Cf2py intent(in)  YIN
278Cf2py depend(NIN) YIN
279Cf2py intent(in) NOUT
280Cf2py intent(in)   XOUT
281Cf2py depend(NOUT) XOUT
282Cf2py intent(out)  YOUT
283Cf2py depend(NOUT) YOUT
284
285      INTEGER NIN,NOUT
286      REAL XIN(NIN),YIN(NIN)
287      REAL XOUT(NOUT),YOUT(NOUT)
288      INTEGER IERROR
289      REAL X,F
290      INTEGER IIN,I
291
292      IERROR = 1
293      IIN = 1
294      X = XOUT(1)
295      DO I=1,NOUT
296         IF (X .GT. XOUT(I)) STOP ! test if Xout not sorted
297         X = XOUT(I)
298         IF (X .LT. XIN(1) .OR. X .GT. XIN(NIN) ) THEN
299            YOUT(I) = 0.0
300         ELSE
301            DO WHILE (X .GT.  XIN(IIN+1))
302               IF (XIN(IIN) .GT. XIN(IIN+1)) STOP ! test if Xin not sorted
303               IIN = IIN + 1
304             ENDDO
305             F = (X - XIN(IIN)) / (XIN(IIN+1) - XIN(IIN))
306             YOUT(I) = (1.-F)*YIN(IIN) + F*YIN(IIN+1)
307          ENDIF
308          !write (*,*) xout(i),iin,f,yout(i)
309      END DO
310      IERROR = 0
311      RETURN
312      END
Note: See TracBrowser for help on using the repository browser.