source: trunk/fsource/pypowder.for @ 415

Last change on this file since 415 was 415, checked in by toby, 11 years ago

implement refinement of constrained parameters; fix minor bugs; add fast interpolate

File size: 5.0 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      FW = (2.355*SQRT(SIG)+GAM)/100.0
18      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
19      FMAX = 15.0*FW
20      DO I=0,NPTS-1
21        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
22     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART)
23      END DO
24      RETURN
25      END
26
27      SUBROUTINE PYDPSVFCJ(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC,
28     1  DPRDT,SIGPART,GAMPART,SLPART)
29C DTT in degrees
30C TTHETA in degrees
31C SPH is S/L + H/L
32C RETURNS FUNCTION & DERIVATIVES
33Cf2py intent(in) NPTS
34Cf2py intent(in) DTT
35cf2py depend(NPTS) DTT
36Cf2py intent(in) TTHETA
37Cf2py intent(in) SIG
38Cf2py intent(in) GAM
39Cf2py intent(in) SPH
40Cf2py intent(out) PRFUNC
41Cf2py depend(NPTS) PRFUNC
42Cf2py intent(out) DPRDT
43Cf2py depend(NPTS) DPRDT
44Cf2py intent(out) SIGPART
45Cf2py depend(NPTS) SIGPART
46Cf2py intent(out) GAMPART
47Cf2py depend(NPTS) GAMPART
48Cf2py intent(out) SLPART
49Cf2py depend(NPTS) SLPART
50
51      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
52     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
53      FW = (2.355*SQRT(SIG)+GAM)/100.0
54      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
55      FMAX = 15.0*FW
56      DO I=0,NPTS-1
57        CALL PSVFCJ(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH,
58     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SLPART(I))
59        DPRDT(I) = DPRDT(I)*100.
60      END DO
61      RETURN
62      END
63
64      SUBROUTINE PYPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SPH,PRFUNC)
65C DTT in degrees
66C TTHETA in degrees
67C SPH is S/L + H/L
68C RETURNS FUNCTION ONLY
69Cf2py intent(in) NPTS
70Cf2py intent(in) DTT
71cf2py depend(NPTS) DTT
72Cf2py intent(in) TTHETA
73Cf2py intent(in) SIG
74Cf2py intent(in) GAM
75Cf2py intent(in) SPH
76Cf2py intent(out) PRFUNC
77Cf2py depend(NPTS) PRFUNC
78
79      REAL*4 DTT(0:NPTS-1),PRFUNC(0:NPTS-1)
80      FW = (2.355*SQRT(SIG)+GAM)/100.0
81      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
82      FMAX = 15.0*FW
83      DO I=0,NPTS-1
84        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SPH/2.0,SPH/2.0,
85     1    PRFUNC(I),DPRDT,SIGPART,GAMPART,SLPART,HLPART)
86      END DO
87      RETURN
88      END
89
90      SUBROUTINE PYDPSVFCJO(NPTS,DTT,TTHETA,SIG,GAM,SHL,PRFUNC,
91     1  DPRDT,SIGPART,GAMPART,SLPART)
92C DTT in degrees
93C TTHETA in degrees
94C SPH is S/L + H/L
95C RETURNS FUNCTION & DERIVATIVES
96Cf2py intent(in) NPTS
97Cf2py intent(in) DTT
98cf2py depend(NPTS) DTT
99Cf2py intent(in) TTHETA
100Cf2py intent(in) SIG
101Cf2py intent(in) GAM
102Cf2py intent(in) SHL
103Cf2py intent(out) PRFUNC
104Cf2py depend(NPTS) PRFUNC
105Cf2py intent(out) DPRDT
106Cf2py depend(NPTS) DPRDT
107Cf2py intent(out) SIGPART
108Cf2py depend(NPTS) SIGPART
109Cf2py intent(out) GAMPART
110Cf2py depend(NPTS) GAMPART
111Cf2py intent(out) SLPART
112Cf2py depend(NPTS) SLPART
113
114      REAL*4 DTT(0:NPTS-1),DPRDT(0:NPTS-1),SIGPART(0:NPTS-1),
115     1  GAMPART(0:NPTS-1),SLPART(0:NPTS-1),PRFUNC(0:NPTS-1)
116      FW = (2.355*SQRT(SIG)+GAM)/100.0
117      FMIN = 10.0*(-FW-SPH*COSD(TTHETA))
118      FMAX = 15.0*FW
119      DO I=0,NPTS-1
120        CALL PSVFCJO(DTT(I)*100.,TTHETA*100.,SIG,GAM,SHL/2.,SHL/2.,
121     1    PRFUNC(I),DPRDT(I),SIGPART(I),GAMPART(I),SPART,HPART)
122          SLPART(I) = SPART
123        DPRDT(I) = DPRDT(I)*100.
124      END DO
125      RETURN
126      END
127
128C Fortran (fast) linear interpolation -- B.H. Toby 9/2011
129      SUBROUTINE PYFINTERP(NIN,XIN,YIN,NOUT,XOUT,YOUT)
130C XIN(1:NIN) and YIN(1:NIN) are arrays of (x,y) points to be interpolated
131C     Values must be sorted increasing in XIN
132C XOUT(1:NOUT) is an array of x values, must also be sorted increasing in x
133C    XOUT may contain values smaller than XIN(1) or larger than XIN(NIN)
134C RETURNS interpolated y values corresponding to XOUT. Values outside the
135C   range of XIN are set to zero.
136C Needs a way to signal an error if XIN or XOUT is not sorted -- for now stops
137Cf2py intent(in) NIN
138Cf2py intent(in)  XIN
139cf2py depend(NIN) XIN
140Cf2py intent(in)  YIN
141cf2py depend(NIN) YIN
142Cf2py intent(in) NOUT
143Cf2py intent(in)   XOUT
144cf2py depend(NOUT) XOUT
145Cf2py intent(out)  YOUT
146cf2py depend(NOUT) YOUT
147
148      REAL XIN(NIN),YIN(NIN)
149      REAL XOUT(NOUT),YOUT(NOUT)
150      INTEGER IERROR
151      REAL X,F
152      INTEGER IIN,I
153
154      IERROR = 1
155      IIN = 1
156      X = XOUT(1)
157      DO I=1,NOUT
158         IF (X .GT. XOUT(I)) STOP ! test if Xout not sorted
159         X = XOUT(I)
160         IF (X .LT. XIN(1) .OR. X .GT. XIN(NIN) ) THEN
161            YOUT(I) = 0.0
162         ELSE
163            DO WHILE (X .GT.  XIN(IIN+1))
164               IF (XIN(IIN) .GT. XIN(IIN+1)) STOP ! test if Xin not sorted
165               IIN = IIN + 1
166             ENDDO
167             F = (X - XIN(IIN)) / (XIN(IIN+1) - XIN(IIN))
168             YOUT(I) = (1.-F)*YIN(IIN) + F*YIN(IIN+1)
169          ENDIF
170          !write (*,*) xout(i),iin,f,yout(i)
171      END DO
172      IERROR = 0
173      RETURN
174      END
Note: See TracBrowser for help on using the repository browser.