source: trunk/fsource/pypowder.for @ 1051

Last change on this file since 1051 was 1051, checked in by vondreele, 10 years ago

remove all mcsa fortran code - not needed
replace 32-bit Windows libraries
complete rework of MC/SA structure factor calcs - no loops!
all numpy/python significant speedup >6x over fortran & >100x over
reflection looped numpy/python
File size: 6.5 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     
184C Fortran (fast) linear interpolation -- B.H. Toby 9/2011
185      SUBROUTINE PYFINTERP(NIN,XIN,YIN,NOUT,XOUT,YOUT)
186C XIN(1:NIN) and YIN(1:NIN) are arrays of (x,y) points to be interpolated
187C     Values must be sorted increasing in XIN
188C XOUT(1:NOUT) is an array of x values, must also be sorted increasing in x
189C    XOUT may contain values smaller than XIN(1) or larger than XIN(NIN)
190C RETURNS interpolated y values corresponding to XOUT. Values outside the
191C   range of XIN are set to zero.
192C Needs a way to signal an error if XIN or XOUT is not sorted -- for now stops
193Cf2py intent(in) NIN
194Cf2py intent(in)  XIN
195cf2py depend(NIN) XIN
196Cf2py intent(in)  YIN
197cf2py depend(NIN) YIN
198Cf2py intent(in) NOUT
199Cf2py intent(in)   XOUT
200cf2py depend(NOUT) XOUT
201Cf2py intent(out)  YOUT
202cf2py depend(NOUT) YOUT
203
204      INTEGER NIN,NOUT
205      REAL XIN(NIN),YIN(NIN)
206      REAL XOUT(NOUT),YOUT(NOUT)
207      INTEGER IERROR
208      REAL X,F
209      INTEGER IIN,I
210
211      IERROR = 1
212      IIN = 1
213      X = XOUT(1)
214      DO I=1,NOUT
215         IF (X .GT. XOUT(I)) STOP ! test if Xout not sorted
216         X = XOUT(I)
217         IF (X .LT. XIN(1) .OR. X .GT. XIN(NIN) ) THEN
218            YOUT(I) = 0.0
219         ELSE
220            DO WHILE (X .GT.  XIN(IIN+1))
221               IF (XIN(IIN) .GT. XIN(IIN+1)) STOP ! test if Xin not sorted
222               IIN = IIN + 1
223             ENDDO
224             F = (X - XIN(IIN)) / (XIN(IIN+1) - XIN(IIN))
225             YOUT(I) = (1.-F)*YIN(IIN) + F*YIN(IIN+1)
226          ENDIF
227          !write (*,*) xout(i),iin,f,yout(i)
228      END DO
229      IERROR = 0
230      RETURN
231      END
Note: See TracBrowser for help on using the repository browser.