source: trunk/fsource/pypowder.for @ 969

Last change on this file since 969 was 969, checked in by vondreele, 9 years ago
File size: 7.8 KB
Line 
1      SUBROUTINE PYPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
2C DTT in degrees
3C TTHETA in degrees
4C SPH is S/L + H/L
5C RETURNS FUNCTION ONLY
6Cf2py intent(in) NPTS
7Cf2py intent(in) DTT
8cf2py depend(NPTS) DTT
9Cf2py intent(in) TTHETA
10Cf2py intent(in) SIG
11Cf2py intent(in) GAM
12Cf2py intent(in) SPH
13Cf2py intent(out) PRFUNC
14Cf2py depend(NPTS) PRFUNC
15
16      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
17      REAL*4 TTHETA,SIG,GAM,SPH
18      INTEGER*4 NPTS,I
19      DO I=0,NPTS-1
20        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
21     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART)
22      END DO
23      RETURN
24      END
25
26      SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC,
27     1  DPRDT,SIGPART,GAMPART,SLPART)
28C DTT in degrees
29C TTHETA in degrees
30C SPH is S/L + H/L
31C RETURNS FUNCTION & DERIVATIVES
32Cf2py intent(in) NPTS
33Cf2py intent(in) DTT
34cf2py depend(NPTS) DTT
35Cf2py intent(in) TTHETA
36Cf2py intent(in) SIG
37Cf2py intent(in) GAM
38Cf2py intent(in) SPH
39Cf2py intent(out) PRFUNC
40Cf2py depend(NPTS) PRFUNC
41Cf2py intent(out) DPRDT
42Cf2py depend(NPTS) DPRDT
43Cf2py intent(out) SIGPART
44Cf2py depend(NPTS) SIGPART
45Cf2py intent(out) GAMPART
46Cf2py depend(NPTS) GAMPART
47Cf2py intent(out) SLPART
48Cf2py depend(NPTS) SLPART
49
50      INTEGER*4 NPTS
51      REAL*4 TTHETA,SIG,GAM,SPH
52      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
53     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
54      DO I=0,NPTS-1
55        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
56     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I))
57        DPRDT(I) = DPRDT(I)*100.
58      END DO
59      RETURN
60      END
61
62      SUBROUTINE PYPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
63C DTT in degrees
64C TTHETA in degrees
65C SPH is S/L + H/L
66C RETURNS FUNCTION ONLY
67Cf2py intent(in) NPTS
68Cf2py intent(in) DTT
69cf2py depend(NPTS) DTT
70Cf2py intent(in) TTHETA
71Cf2py intent(in) SIG
72Cf2py intent(in) GAM
73Cf2py intent(in) SPH
74Cf2py intent(out) PRFUNC
75Cf2py depend(NPTS) PRFUNC
76
77      INTEGER*4 NPTS
78      REAL*4 TTHETA,SIG,GAM,SPH
79      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
80      DO I=0,NPTS-1
81        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
82     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
83      END DO
84      RETURN
85      END
86
87      SUBROUTINE PYDPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SHL,PRFUNC,
88     1  DPRDT,SIGPART,GAMPART,SLPART)
89C DTT in degrees
90C TTHETA in degrees
91C SPH is S/L + H/L
92C RETURNS FUNCTION & DERIVATIVES
93Cf2py intent(in) NPTS
94Cf2py intent(in) DTT
95cf2py depend(NPTS) DTT
96Cf2py intent(in) TTHETA
97Cf2py intent(in) SIG
98Cf2py intent(in) GAM
99Cf2py intent(in) SHL
100Cf2py intent(out) PRFUNC
101Cf2py depend(NPTS) PRFUNC
102Cf2py intent(out) DPRDT
103Cf2py depend(NPTS) DPRDT
104Cf2py intent(out) SIGPART
105Cf2py depend(NPTS) SIGPART
106Cf2py intent(out) GAMPART
107Cf2py depend(NPTS) GAMPART
108Cf2py intent(out) SLPART
109Cf2py depend(NPTS) SLPART
110
111      INTEGER*4 NPTS
112      REAL*4 TTHETA,SIG,GAM,SHL
113      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
114     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
115      DO I=0,NPTS-1
116        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
117     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
118          SLPART(I) = SPART
119        DPRDT(I) = DPRDT(I)*100.
120      END DO
121      RETURN
122      END
123
124      SUBROUTINE PYEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC)
125C DTT in microsec
126C RETURNS FUNCTION
127Cf2py intent(in) NPTS
128Cf2py intent(in) DTT
129cf2py depend(NPTS) DTT
130Cf2py intent(in) ALP
131Cf2py intent(in) BET
132Cf2py intent(in) SIG
133Cf2py intent(in) GAM
134Cf2py intent(out) PRFUNC
135Cf2py depend(NPTS) PRFUNC
136
137      INTEGER*4 NPTS
138      REAL*4 ALP,BET,SIG,GAM,SHL
139      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1),DPRDT,ALPPART,
140     1  BETPART,SIGPART,GAMPART
141      DO I=0,NPTS-1
142        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT,
143     1    ALPPART,BETPART,SIGPART,GAMPART)
144      END DO
145      RETURN
146      END
147
148      SUBROUTINE PYDEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC,
149     1  DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
150C DTT in microsec
151C RETURNS FUNCTION & DERIVATIVES
152Cf2py intent(in) NPTS
153Cf2py intent(in) DTT
154cf2py depend(NPTS) DTT
155Cf2py intent(in) ALP
156Cf2py intent(in) BET
157Cf2py intent(in) SIG
158Cf2py intent(in) GAM
159Cf2py intent(out) PRFUNC
160Cf2py depend(NPTS) PRFUNC
161Cf2py intent(out) DPRDT
162Cf2py depend(NPTS) DPRDT
163Cf2py intent(out) ALPPART
164Cf2py depend(NPTS) ALPPART
165Cf2py intent(out) BETPART
166Cf2py depend(NPTS) BETPART
167Cf2py intent(out) SIGPART
168Cf2py depend(NPTS) SIGPART
169Cf2py intent(out) GAMPART
170Cf2py depend(NPTS) GAMPART
171
172      INTEGER*4 NPTS
173      REAL*4 ALP,BET,SIG,GAM,SHL
174      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),ALPPART(0:NPTS-1),
175     1  BETPART(0:NPTS-1),SIGPART(0:NPTS-1),
176     1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
177      DO I=0,NPTS-1
178        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT(I),
179     1    ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
180      END DO
181      RETURN
182      END
183
184      SUBROUTINE PYMCSASFCALC(INV,NTD,TDATA,MDATA,X,MUL,NFFS,FFS,
185     1  NUNIQ,U,PHI,ICALC)
186Cf2py intent(in) INV
187Cf2py intent(in) NTD
188Cf2py intent(in) TDATA
189cf2py depend(NTD) TDATA
190Cf2py intent(in) MDATA
191cf2py depend(NTD) MDATA
192Cf2py intent(in) X
193cf2py depend(NTD) X
194Cf2py intent(in) MUL
195Cf2py intent(in) NFFS
196Cf2py intent(in) FFS
197cf2py depend(NFFS) FFS
198Cf2py intent(in) NUNIQ
199Cf2py intent(in) U
200cf2py depend(NUNIQ) U
201Cf2py intent(in) PHI
202cf2py depend(NUNIQ) PHI
203Cf2py intent(out) ICALC
204
205      LOGICAL*4 INV
206      INTEGER*4 NTD,MUL,NFFS,NUNIQ,I,J,K,TDATA(0:NTD-1)
207      REAL*4 X(0:3*NTD-1),U(0:3*NUNIQ-1)
208      REAL*4 MDATA(0:NTD-1),FFS(0:NFFS-1)
209      REAL*4 ICALC,PHI(0:NUNIQ-1)
210      REAL*4 PHASE(0:NTD,0:NUNIQ),FF(0:NTD-1),FAS,FBS,TWOPI
211
212      TWOPI = 8.0*ATAN(1.0)
213      DO I=0,NTD-1
214        FF(I) = FFS(TDATA(I))*MDATA(I)/NUNIQ
215        DO J=0,NUNIQ-1
216          PHASE(I,J) = 0.
217          DO K=0,2
218            PHASE(I,J) = PHASE(I,J)+U(3*J+K)*X(3*I+K)
219          END DO
220          PHASE(I,J) = PHASE(I,J)+PHI(J)
221        END DO
222      END DO
223      FAS = 0.
224      FBS = 0.
225      DO I=0,NTD-1
226        DO J=0,NUNIQ-1
227          FAS = FAS+FF(I)*COS(TWOPI*PHASE(I,J))
228          IF ( .NOT. INV ) FBS = FBS+FF(I)*SIN(TWOPI*PHASE(I,J))
229        END DO
230      END DO
231      ICALC = (FAS**2+FBS**2)*MUL
232      RETURN
233      END
234
235C Fortran (fast) linear interpolation -- B.H. Toby 9/2011
236      SUBROUTINE PYFINTERP(NIN,XIN,YIN,NOUT,XOUT,YOUT)
237C XIN(1:NIN) and YIN(1:NIN) are arrays of (x,y) points to be interpolated
238C     Values must be sorted increasing in XIN
239C XOUT(1:NOUT) is an array of x values, must also be sorted increasing in x
240C    XOUT may contain values smaller than XIN(1) or larger than XIN(NIN)
241C RETURNS interpolated y values corresponding to XOUT. Values outside the
242C   range of XIN are set to zero.
243C Needs a way to signal an error if XIN or XOUT is not sorted -- for now stops
244Cf2py intent(in) NIN
245Cf2py intent(in)  XIN
246cf2py depend(NIN) XIN
247Cf2py intent(in)  YIN
248cf2py depend(NIN) YIN
249Cf2py intent(in) NOUT
250Cf2py intent(in)   XOUT
251cf2py depend(NOUT) XOUT
252Cf2py intent(out)  YOUT
253cf2py depend(NOUT) YOUT
254
255      INTEGER NIN,NOUT
256      REAL XIN(NIN),YIN(NIN)
257      REAL XOUT(NOUT),YOUT(NOUT)
258      INTEGER IERROR
259      REAL X,F
260      INTEGER IIN,I
261
262      IERROR = 1
263      IIN = 1
264      X = XOUT(1)
265      DO I=1,NOUT
266         IF (X .GT. XOUT(I)) STOP ! test if Xout not sorted
267         X = XOUT(I)
268         IF (X .LT. XIN(1) .OR. X .GT. XIN(NIN) ) THEN
269            YOUT(I) = 0.0
270         ELSE
271            DO WHILE (X .GT.  XIN(IIN+1))
272               IF (XIN(IIN) .GT. XIN(IIN+1)) STOP ! test if Xin not sorted
273               IIN = IIN + 1
274             ENDDO
275             F = (X - XIN(IIN)) / (XIN(IIN+1) - XIN(IIN))
276             YOUT(I) = (1.-F)*YIN(IIN) + F*YIN(IIN+1)
277          ENDIF
278          !write (*,*) xout(i),iin,f,yout(i)
279      END DO
280      IERROR = 0
281      RETURN
282      END
Note: See TracBrowser for help on using the repository browser.