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