source: trunk/fsource/pypowder.for @ 1274

Last change on this file since 1274 was 1274, checked in by vondreele, 8 years ago

implement Bragg peaks in SASD - required new fortran routines
implement unified (Guinier + Porod) and Porod models in SASD
implement monodisperse models in SASD
correct bin width issue in lognormal, etc. fitting

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 PYPSVOIGT(NPTS,DTT,SIG,GAM,PRFUNC)
63C DTT in degrees
64C RETURNS FUNCTION ONLY
65Cf2py intent(in) NPTS
66Cf2py intent(in) DTT
67cf2py depend(NPTS) DTT
68Cf2py intent(in) SIG
69Cf2py intent(in) GAM
70Cf2py intent(out) PRFUNC
71Cf2py depend(NPTS) PRFUNC
72
73      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
74      REAL*4 SIG,GAM
75      INTEGER*4 NPTS,I
76      DO I=0,NPTS-1
77        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
78     1    PRFUNC(I),DPRDT,SIGPART,GAMPART)
79      END DO
80      RETURN
81      END
82
83      SUBROUTINE PYDPSVOIGT(NPTS,DTT,SIG,GAM,PRFUNC,
84     1  DPRDT,SIGPART,GAMPART)
85C DTT in degrees
86C RETURNS FUNCTION & DERIVATIVES
87Cf2py intent(in) NPTS
88Cf2py intent(in) DTT
89cf2py depend(NPTS) DTT
90Cf2py intent(in) SIG
91Cf2py intent(in) GAM
92Cf2py intent(out) PRFUNC
93Cf2py depend(NPTS) PRFUNC
94Cf2py intent(out) DPRDT
95Cf2py depend(NPTS) DPRDT
96Cf2py intent(out) SIGPART
97Cf2py depend(NPTS) SIGPART
98Cf2py intent(out) GAMPART
99Cf2py depend(NPTS) GAMPART
100
101      INTEGER*4 NPTS
102      REAL*4 SIG,GAM
103      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
104     1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
105      DO I=0,NPTS-1
106        CALL PSVOIGT(DTT(I)*100.,SIG,GAM,
107     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I))
108        DPRDT(I) = DPRDT(I)*100.
109      END DO
110      RETURN
111      END
112
113      SUBROUTINE PYPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
114C DTT in degrees
115C TTHETA in degrees
116C SPH is S/L + H/L
117C RETURNS FUNCTION ONLY
118Cf2py intent(in) NPTS
119Cf2py intent(in) DTT
120cf2py depend(NPTS) DTT
121Cf2py intent(in) TTHETA
122Cf2py intent(in) SIG
123Cf2py intent(in) GAM
124Cf2py intent(in) SPH
125Cf2py intent(out) PRFUNC
126Cf2py depend(NPTS) PRFUNC
127
128      INTEGER*4 NPTS
129      REAL*4 TTHETA,SIG,GAM,SPH
130      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
131      DO I=0,NPTS-1
132        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
133     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
134      END DO
135      RETURN
136      END
137
138      SUBROUTINE PYDPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SHL,PRFUNC,
139     1  DPRDT,SIGPART,GAMPART,SLPART)
140C DTT in degrees
141C TTHETA in degrees
142C SPH is S/L + H/L
143C RETURNS FUNCTION & DERIVATIVES
144Cf2py intent(in) NPTS
145Cf2py intent(in) DTT
146cf2py depend(NPTS) DTT
147Cf2py intent(in) TTHETA
148Cf2py intent(in) SIG
149Cf2py intent(in) GAM
150Cf2py intent(in) SHL
151Cf2py intent(out) PRFUNC
152Cf2py depend(NPTS) PRFUNC
153Cf2py intent(out) DPRDT
154Cf2py depend(NPTS) DPRDT
155Cf2py intent(out) SIGPART
156Cf2py depend(NPTS) SIGPART
157Cf2py intent(out) GAMPART
158Cf2py depend(NPTS) GAMPART
159Cf2py intent(out) SLPART
160Cf2py depend(NPTS) SLPART
161
162      INTEGER*4 NPTS
163      REAL*4 TTHETA,SIG,GAM,SHL
164      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
165     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
166      DO I=0,NPTS-1
167        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
168     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
169          SLPART(I) = SPART
170        DPRDT(I) = DPRDT(I)*100.
171      END DO
172      RETURN
173      END
174
175      SUBROUTINE PYEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC)
176C DTT in microsec
177C RETURNS FUNCTION
178Cf2py intent(in) NPTS
179Cf2py intent(in) DTT
180cf2py depend(NPTS) DTT
181Cf2py intent(in) ALP
182Cf2py intent(in) BET
183Cf2py intent(in) SIG
184Cf2py intent(in) GAM
185Cf2py intent(out) PRFUNC
186Cf2py depend(NPTS) PRFUNC
187
188      INTEGER*4 NPTS
189      REAL*4 ALP,BET,SIG,GAM,SHL
190      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1),DPRDT,ALPPART,
191     1  BETPART,SIGPART,GAMPART
192      DO I=0,NPTS-1
193        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT,
194     1    ALPPART,BETPART,SIGPART,GAMPART)
195      END DO
196      RETURN
197      END
198
199      SUBROUTINE PYDEPSVOIGT(NPTS,DTT,ALP,BET,SIG,GAM,PRFUNC,
200     1  DPRDT,ALPPART,BETPART,SIGPART,GAMPART)
201C DTT in microsec
202C RETURNS FUNCTION & DERIVATIVES
203Cf2py intent(in) NPTS
204Cf2py intent(in) DTT
205cf2py depend(NPTS) DTT
206Cf2py intent(in) ALP
207Cf2py intent(in) BET
208Cf2py intent(in) SIG
209Cf2py intent(in) GAM
210Cf2py intent(out) PRFUNC
211Cf2py depend(NPTS) PRFUNC
212Cf2py intent(out) DPRDT
213Cf2py depend(NPTS) DPRDT
214Cf2py intent(out) ALPPART
215Cf2py depend(NPTS) ALPPART
216Cf2py intent(out) BETPART
217Cf2py depend(NPTS) BETPART
218Cf2py intent(out) SIGPART
219Cf2py depend(NPTS) SIGPART
220Cf2py intent(out) GAMPART
221Cf2py depend(NPTS) GAMPART
222
223      INTEGER*4 NPTS
224      REAL*4 ALP,BET,SIG,GAM,SHL
225      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),ALPPART(0:NPTS-1),
226     1  BETPART(0:NPTS-1),SIGPART(0:NPTS-1),
227     1  GAMPART(0:NPTS-1),PRFUNC(0:NPTS-1)
228      DO I=0,NPTS-1
229        CALL EPSVOIGT(DTT(I),ALP,BET,SIG,GAM,PRFUNC(I),DPRDT(I),
230     1    ALPPART(I),BETPART(I),SIGPART(I),GAMPART(I))
231      END DO
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.