V3FIT
furpl.f
1  SUBROUTINE furplm(p,grsq,x,y,v,xp,yp,gsq,np,nx,ny,ntx,nw,nh,is,js)
2 c ***********************************************************************
3 c given an x-y grid with a variable p defined on the grid this routine
4 c finds the level p curves of value v. that is, it finds np points
5 c (xp(k),yp(k)) such that p(xp(k),yp(k)) is CLOSE to v.
6 c id is an INTEGER between 1 and 12 or 15 which identifies the option
7 c USEd to find the last point.
8 c IF in=1 the curve is moving up and IF in=0 the curve is moving DOwn.
9 c IF il=1 the curve is moving left and IF il=0 the curve is moving right.
10 c IF the last poINT was a grid poINT id=15 and igp is an INTEGER between
11 c 1 and 12 which identifies the grid poINT option USEd to find the last
12 c point.
13 c IF ig=1 the grid poINT was entered from the left and IF ig=0 the grid
14 c poINT was entered from the right.
15 c iq is an indicator USEd in that part of the PROGRAM which searches the
16 c grid squares surrounding a grid point.
17 c IF ntx is less than zero search for the first poINT starts from the
18 c top of the grid and goes DOwn.
19 c IF np eq -1 the search starts in the center of the grid
20 c IF ntx gt 0 the search starts at the bottom of the grid
21 c IF you would like further details see joanne helton.
22 c ***********************************************************************
23  USE precision
24  IMPLICIT NONE
25  INTEGER :: js, is, ntx, ny, nx, np, nw, nh, nt, nx2, ny2
26  INTEGER :: jst, jdec, kprt, j, i, k, in, iq, id, i1, jy,
27  . ix, igp, l, iyy, ixx, ig, il
28  REAL(rprec) :: xp(*),yp(*)
29  REAL(rprec) :: grsq(nw,nh),gsq(*)
30  REAL(rprec) :: p(nw,*),x(*),y(*)
31  REAL(rprec) :: a, b, c, f, yq, v, dx, dy, dp, d, e, ff, ws,
32  . xl, yl, x1, y1
33  a=0.
34  b=0.
35  c=0.
36  f=0.
37  yq=0.
38  nt=ntx
39  jst=js
40  jdec=1
41  IF(np.ne.-1) go to 1999
42  jst=ny/2+1
43  jdec=1
44  nt=iabs(nt)
45  1999 IF(nt.gt.0) go to 2004
46  jst=ny-1
47  jdec=-1
48  nt=-nt
49  2004 CONTINUE
50  nx2=nx/2
51  ny2=ny/2
52  dx=x(2)-x(1)
53  dy=y(2)-y(1)
54  dp=abs(p(nx2,ny2))
55  d=10.**(-1*8)*dp
56  e=10.**(-1*4)*dx*3.0
57  ff=10.**(-1*6)*dx
58  kprt=0
59 2001 CONTINUE
60  j=jst
61  i=is
62  k=0
63  in=0
64  iq=0
65 c ---------------------
66 c we find a first point
67 c ---------------------
68  13 IF((v-p(i,j))*(p(i+1,j)-v).gt.0.) go to 30
69  IF(abs(p(i,j)-v).ge.d) go to 32
70  in=1
71  ig=1
72  go to 15
73  32 i=i+1
74  IF(i.eq.nx) go to 25
75  go to 13
76  25 j=j+jdec
77  IF(jdec.lt.0) go to 2005
78  IF(j.gt.ny) go to 10
79  i=1
80  go to 13
81  2005 IF(j.lt.1) go to 10
82  i=1
83  go to 13
84  30 IF(i.gt.(nx-2)) go to 74
85  IF(i.eq.1) go to 88
86  xp(1)=x(i)
87  CALL fit(2,x(i-1),x(i),x(i+1),x(i+2),p(i-1,j),p(i,j),p(i+1,j),
88  1 p(i+2,j),xp(1),v,yq)
89  ws=xp(1)
90  CALL cubic(x(i-1),x(i),x(i+1),x(i+2),ws,a,b,c,f)
91  gsq(1)=a*grsq(i-1,j)+b*grsq(i,j)+c*grsq(i+1,j)+f*grsq(i+2,j)
92  go to 75
93  88 xp(1)=x(2)
94  CALL fit(2,x(1),x(2),x(3),x(4),p(1,j),p(2,j),p(3,j),p(4,j),
95  1 xp(1),v,yq)
96  ws=xp(1)
97  CALL cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
98  gsq(1)=a*grsq(1,j)+b*grsq(2,j)+c*grsq(3,j)+f*grsq(4,j)
99  go to 75
100  74 xp(1)=x(nx-2)
101  CALL fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,j),p(nx-2,j),
102  1 p(nx-1,j),p(nx,j),xp(1),v,yq)
103  ws=xp(1)
104  CALL cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
105  gsq(1)=a*grsq(nx-3,j)+b*grsq(nx-2,j)+c*grsq(nx-1,j)+f*grsq(nx,j)
106  75 yp(1)=y(j)
107  in=1
108  IF(xp(1).le.x(i)) xp(1)=x(i)+ff
109  IF(xp(1).ge.x(i+1)) xp(1)=x(i+1)-ff
110  id=20
111  il=1
112  IF(jdec.lt.0) il=0
113  k=1
114  IF(kprt.eq.1) WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
115  x ,iq,il
116  102 FORMAT(i5,3e13.6,10i5)
117  go to 5
118 c ---------------------------------------------------------------
119 c in the block below we determine WHERE the next poINT is located
120 c ---------------------------------------------------------------
121 c ---------------------------------------------------------------------
122 c the curve is entering the grid square from the right or from the left
123 c IF il equals 1 the curve is moving left
124 c ---------------------------------------------------------------------
125  4 IF(il.eq.1) go to 21
126 c --------------------------------------------------
127 c the following options are for a curve moving right
128 c --------------------------------------------------
129  IF(i+1.gt.nx) go to 10
130 c --------------------------
131 c first we check grid points
132 c --------------------------
133  IF(abs(v-p(i+1,j+1)).ge.d) go to 45
134  ig=1
135  in=1
136  i=i+1
137  j=j+1
138  igp=11
139  go to 15
140  45 IF(abs(v-p(i+1,j)).ge.d) go to 33
141  ig=1
142  in=0
143  i=i+1
144  igp=12
145  go to 15
146 c --------------------------------
147 c now we check between grid points
148 c --------------------------------
149  33 IF((v-p(i,j))*(p(i+1,j)-v).le.0.) go to 50
150  ix=i
151  jy=j
152  id=1
153  go to 1
154  50 IF((v-p(i,j+1))*(p(i+1,j+1)-v).le.0) go to 51
155  ix=i
156  jy=j+1
157  id=2
158  go to 1
159  51 IF((v-p(i+1,j))*(p(i+1,j+1)-v).le.0.) go to 52
160  ix=i+1
161  jy=j
162  id=3
163  go to 7
164  52 CONTINUE
165 c -------------------------------------------------
166 c the following options are for a curve moving left
167 c -------------------------------------------------
168  IF(i-1.lt.0) go to 10
169 c --------------------------
170 c first we check grid points
171 c --------------------------
172  21 IF(abs(v-p(i-1,j)).ge.d) go to 46
173  in=0
174  ig=0
175  i=i-1
176  igp=1
177  go to 15
178  46 IF(abs(v-p(i-1,j+1)).ge.d) go to 34
179  ig=0
180  in=1
181  i=i-1
182  j=j+1
183  igp=2
184  go to 15
185 c --------------------------------
186 c now we check between grid points
187 c --------------------------------
188  34 IF((v-p(i,j+1))*(p(i-1,j+1)-v).le.0.) go to 53
189  ix=i-1
190  jy=j+1
191  id=4
192  go to 1
193  53 IF((v-p(i-1,j+1))*(p(i-1,j)-v).le.0) go to 54
194  ix=i-1
195  jy=j
196  id=5
197  go to 7
198  54 IF((v-p(i,j))*(p(i-1,j)-v).le.0.) go to 55
199  ix=i-1
200  jy=j
201  id=6
202  go to 1
203  55 CONTINUE
204 c ---------------------------------------------------------------------
205 c the curve is entering the grid square from the top or from the bottom
206 c or from a grid point
207 c IF ig equals 1 the grid poINT was entered from the left
208 c IF in equals 1 the curve is moving up
209 c ---------------------------------------------------------------------
210  5 IF(id.eq.15.and.il.eq.1) i=i-1
211  72 iq=0
212  62 IF(in.eq.1) go to 12
213 c -------------------------------------------------
214 c the following options are for a curve moving DOwn
215 c -------------------------------------------------
216  IF(j-1.lt.1) go to 10
217  11 IF(abs(v-p(i+1,j-1)).ge.d) go to 42
218  ig=1
219  in=0
220  i=i+1
221  j=j-1
222  igp=3
223  go to 15
224  42 IF(abs(v-p(i,j-1)).ge.d) go to 17
225  ig=0
226  in=0
227  j=j-1
228  igp=5
229  go to 15
230 c --------------------------------
231 c now we check between grid points
232 c --------------------------------
233  17 IF((v-p(i+1,j-1))*(p(i+1,j)-v).le.0.) go to 56
234  IF(abs(v-p(i+1,j)).le.d) go to 56
235  IF(id.eq.15.and.igp.eq.7) go to 56
236  ix=i+1
237  jy=j-1
238  id=7
239  go to 7
240  56 IF((v-p(i,j-1))*(p(i+1,j-1)-v).le.0.) go to 57
241  ix=i
242  jy=j-1
243  id=8
244  go to 1
245  57 IF((v-p(i,j))*(p(i,j-1)-v).le.0.) go to 44
246  IF(abs(v-p(i,j)).le.d) go to 44
247  ix=i
248  jy=j-1
249  id=9
250  go to 7
251  44 IF(iq.ne.1) go to 48
252  in=1
253  IF(il.eq.1) i=i+1
254  IF(il.eq.0) i=i-1
255  go to 72
256  48 iq=1
257  in=1
258  go to 62
259 c -----------------------------------------------
260 c the following options are for a curve moving up
261 c -----------------------------------------------
262  12 IF(j+1.gt.ny) go to 10
263  IF(i.eq.nx) go to 10
264 c --------------------------
265 c first we check grid points
266 c --------------------------
267  IF(abs(v-p(i+1,j+1)).ge.d) go to 40
268  ig=1
269  in=1
270  i=i+1
271  j=j+1
272  igp=7
273  go to 15
274  40 IF(abs(v-p(i,j+1)).ge.d) go to 61
275  ig=0
276  in=1
277  j=j+1
278  igp=9
279  go to 15
280 c --------------------------------
281 c now we check between grid points
282 c --------------------------------
283  61 IF((v-p(i,j+1))*(p(i+1,j+1)-v).le.0) go to 59
284  ix=i
285  jy=j+1
286  id=10
287  go to 1
288  59 IF((v-p(i+1,j))*(p(i+1,j+1)-v).le.0.) go to 60
289  IF(abs(v-p(i+1,j)).le.d) go to 60
290  ix=i+1
291  jy=j
292  id=11
293  go to 7
294  60 IF((v-p(i,j+1))*(p(i,j)-v).le.0.) go to 41
295  IF(abs(v-p(i,j)).le.d) go to 41
296  ix=i
297  jy=j
298  id=12
299  go to 7
300  41 IF(iq.ne.1) go to 49
301  in=0
302  IF(il.eq.1) i=i+1
303  IF(il.eq.0) i=i-1
304  go to 72
305  49 iq=1
306  in=0
307  go to 62
308 c ---------------------------------------------------------------------
309 c in the block below the x and y values of the poINT found are computed
310 c ---------------------------------------------------------------------
311 c and stored and control is RETURNed to the appropriate poINT in the
312 c block above
313  1 k=k+1
314  IF(k.eq.nt) WRITE(66,106)
315  106 FORMAT(/,"***nt points found in furplm***",/)
316  IF(k.eq.nt) go to 2000
317  IF(ix.eq.1) go to 2
318  IF(ix.gt.(nx-2)) go to 3
319  xp(k)=x(ix)
320  CALL fit(2,x(ix-1),x(ix),x(ix+1),x(ix+2),p(ix-1,jy),p(ix,jy),
321  1 p(ix+1,jy),p(ix+2,jy),xp(k),v,yq)
322  ws=xp(k)
323  CALL cubic(x(ix-1),x(ix),x(ix+1),x(ix+2),ws,a,b,c,f)
324  2009 CONTINUE
325  gsq(k)=a*grsq(ix-1,jy)+b*grsq(ix,jy)+c*grsq(ix+1,jy)
326  1+f*grsq(ix+2,jy)
327  2010 CONTINUE
328  go to 6
329  2 xp(k)=x(2)
330  CALL fit(2,x(1),x(2),x(3),x(4),p(1,jy),p(2,jy),p(3,jy),p(4,jy),
331  1 xp(k),v,yq)
332  ws=xp(k)
333  CALL cubic(x(1),x(2),x(3),x(4),ws,a,b,c,f)
334  gsq(k)=a*grsq(1,jy)+b*grsq(2,jy)+c*grsq(3,jy)+f*grsq(4,jy)
335  go to 6
336  3 xp(k)=x(nx-2)
337  CALL fit(2,x(nx-3),x(nx-2),x(nx-1),x(nx),p(nx-3,jy),p(nx-2,jy),
338  1 p(nx-1,jy),p(nx,jy),xp(k),v,yq)
339  ws=xp(k)
340  CALL cubic(x(nx-3),x(nx-2),x(nx-1),x(nx),ws,a,b,c,f)
341  gsq(k)=a*grsq(nx-3,jy)+b*grsq(nx-2,jy)+c*grsq(nx-1,jy)
342  1+f*grsq(nx,jy)
343  6 yp(k)=y(jy)
344  IF((id.eq.4).or.(id.eq.6)) i=i-1
345  IF(xp(k).le.x(ix)) xp(k)=x(ix)+ff
346  IF(xp(k).ge.x(ix+1)) xp(k)=x(ix+1)-ff
347  IF(k.le.3) go to 20
348  IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
349  1.and.k.ne.1) go to 10
350  20 IF(jy.eq.1.and.(((id.eq.1).or.(id.eq.6)).or.(id.eq.8))) go to 10
351  IF(jy.eq.ny.and.(((id.eq.2).or.(id.eq.4)).or.(id.eq.10))) go to 10
352  IF(((id.eq.2).or.(id.eq.4)).or.(id.eq.10)) j=j+1
353  IF(id.eq.8) j=j-1
354  in=1
355  IF(yp(k).lt.yp(k-1)) in=0
356  il=1
357  IF(xp(k).gt.xp(k-1)) il=0
358  IF(k > 2) THEN ! avoid reverence to element (0) that does not exist
359  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
360  x 3) WRITE(66,107)
361  107 FORMAT(/,"***poINT same as poINT before last***",/)
362  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
363  1.and.k.ge.3) go to 2000
364  ENDIF
365  IF(kprt.eq.1) WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
366  x ,iq,il
367  go to 5
368  7 k=k+1
369  IF(k.eq.nt) WRITE(66,106)
370  IF(k.eq.nt) go to 2000
371  xp(k)=x(ix)
372  IF(jy.eq.1) go to 8
373  IF(jy.gt.(ny-2)) go to 9
374  yp(k)=y(jy)
375  CALL fit(2,y(jy-1),y(jy),y(jy+1),y(jy+2),p(ix,jy-1),p(ix,jy),
376  1 p(ix,jy+1),p(ix,jy+2),yp(k),v,yq)
377  ws=yp(k)
378  CALL cubic(y(jy-1),y(jy),y(jy+1),y(jy+2),ws,a,b,c,f)
379  gsq(k)=a*grsq(ix,jy-1)+b*grsq(ix,jy)+c*grsq(ix,jy+1)
380  1+f*grsq(ix,jy+2)
381  go to 14
382  8 yp(k)=y(2)
383  CALL fit(2,y(1),y(2),y(3),y(4),p(ix,1),p(ix,2),p(ix,3),p(ix,4),
384  1 yp(k),v,yq)
385  ws=yp(k)
386  CALL cubic(y(1),y(2),y(3),y(4),ws,a,b,c,f)
387  gsq(k)=a*grsq(ix,1)+b*grsq(ix,2)+c*grsq(ix,3)+f*grsq(ix,4)
388  go to 14
389  9 yp(k)=y(ny-2)
390  CALL fit(2,y(ny-3),y(ny-2),y(ny-1),y(ny),p(ix,ny-3),p(ix,ny-2),
391  1 p(ix,ny-1),p(ix,ny),yp(k),v,yq)
392  ws=yp(k)
393  CALL cubic(y(ny-3),y(ny-2),y(ny-1),y(ny),ws,a,b,c,f)
394  gsq(k)=a*grsq(ix,ny-3)+b*grsq(ix,ny-2)+c*grsq(ix,ny-1)
395  1+f*grsq(ix,ny)
396  14 IF(yp(k).le.y(jy)) yp(k)=y(jy)+ff
397  IF(yp(k).ge.y(jy+1)) yp(k)=y(jy+1)-ff
398  IF(k.le.3) go to 19
399  IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
400  1.and.k.ne.1) go to 10
401  19 IF(ix.eq.1.and.(((id.eq.5).or.(id.eq.9)).or.(id.eq.12))) go to 10
402  IF(ix.eq.nx.and.(((id.eq.3).or.(id.eq.7)).or.(id.eq.11))) go to 10
403  IF(id.eq.5) i=i-1
404  IF((id.eq.7).or.(id.eq.9)) j=j-1
405  in=1
406  IF(yp(k).lt.yp(k-1)) in=0
407  il=1
408  IF(xp(k).gt.xp(k-1)) il=0
409  IF(((id.eq.3).or.(id.eq.7.and.il.eq.0)).or.(id.eq.11.and.il.eq.0))
410  1i=i+1
411  IF(kprt.eq.1) WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
412  x ,iq,il
413  IF(k > 2) THEN ! avoid reverence to element (0) that does not exist
414  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
415  x 3) WRITE(66,107)! avoid reverence to element (0) that does not exist
416  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
417  1.and.k.ge.3) go to 2000
418  ENDIF
419  go to 4
420  15 k=k+1
421  IF(k.eq.nt) WRITE(66,106)
422  IF(k.eq.nt) go to 2000
423  xp(k)=x(i)
424  yp(k)=y(j)
425  IF(k > 2) THEN ! avoid reverence to element (0) that does not exist
426  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e).and.k.ge.
427  x 3) WRITE(66,107)
428  IF((abs(xp(k)-xp(k-2)).lt.e.and.abs(yp(k)-yp(k-2)).lt.e)
429  1.and.k.ge.3) go to 2000
430  ENDIF
431  gsq(k)=grsq(i,j)
432  id=15
433  IF(k.le.3) go to 18
434  IF((abs(xp(k)-xp(1)).lt.e.and.abs(yp(k)-yp(1)).lt.e)
435  1.and.k.ne.1) go to 10
436  18 IF(((i.eq.1.or.i.eq.nx).and.(j.eq.1.or.j.eq.ny)).and.k.ne.1)
437  1go to 10
438  in=1
439  IF(yp(k).lt.yp(k-1)) in=0
440  il=1
441  IF(xp(k).gt.xp(k-1)) il=0
442  IF(kprt.eq.1) WRITE(66,102)k,xp(k),yp(k),p(i,j),i,j,ix,jy,ig,in,id
443  x ,iq,il,igp
444  go to 5
445  10 CONTINUE
446  np=k
447  RETURN
448  2000 CONTINUE
449 c test IF looping on separatrix
450  DO 2020 i=1,k
451  x1=xp(i)
452  y1=yp(i)
453  DO 2008 j=i+1,k
454  IF(abs(xp(j)-x1).gt.e) go to 2008
455  IF(abs(yp(j)-y1).le.e) go to 2030
456  2008 CONTINUE
457  2020 CONTINUE
458  go to 2100
459  2030 CONTINUE
460  k=j-i+1
461  DO j=1,k
462  l=i+j-1
463  gsq(j)=gsq(l)
464  xp(j)=xp(l)
465  yp(j)=yp(l)
466  END DO
467  go to 10
468  2100 CONTINUE
469  WRITE(6,*)'error in furpl--fur1'
470  stop
471  WRITE(66,2002) v,(ixx,x(ixx),ixx=1,nx),(iyy,y(iyy),iyy=1,ny)
472 2002 FORMAT(" furplm diagnostic center",/,
473  . " value =",e14.6,/," x and y grids follow",/,
474  .(i5,e14.6,i5,e14.6,i5,e14.6,i5,e14.6))
475  WRITE(66,2003)
476 2003 FORMAT(//," k "," xp yp p(i,j) ",
477  ." i j ix jy ig in id iq il")
478  kprt=1
479  go to 2001
480  END SUBROUTINE furplm
481  SUBROUTINE cubic(p1,p2,p3,p4,p,a,b,c,d)
482  USE precision
483  IMPLICIT NONE
484  REAL(rprec) :: p1, p2, p3, p4, p, a, b, c, d, pp1, pp2, pp3, pp4
485  REAL(rprec) :: zero, p1p2, p2p3, p1p3, p1p4, p2p4, p3p4
486  INTEGER :: k12, k23, k34
487  data zero/1.e-6_dbl/
488 cray lcm (p1),(p2),(p3),(p4)
489  p1p2=p1-p2
490  k12=sign(1.5_dbl,p1p2)
491  IF(abs(p1p2).lt.zero) go to 1
492  p2p3=p2-p3
493  k23=sign(1.5_dbl,p2p3)
494  IF(abs(p2p3).lt.zero) go to 2
495  p1p3=p1-p3
496  IF(abs(p1p3).lt.zero) go to 1
497  p1p4=p1-p4
498  IF (abs(p1p4).lt.zero) go to 3
499  p2p4=p2-p4
500  IF (abs(p2p4).lt.zero) go to 4
501  p3p4=p3-p4
502  k34=sign(1.5_dbl,p3p4)
503  IF(abs(p3p4).lt.zero) go to 4
504  IF (iabs(k12+k23+k34).eq.3) go to 5
505  IF (iabs(k23+k34).eq.2) go to 1
506  IF (iabs(k12+k23).eq.2) go to 4
507 5 pp1=p-p1
508  pp2=p-p2
509  pp3=p-p3
510  pp4=p-p4
511  a=pp2*pp3*pp4/(p1p2*p1p3*p1p4)
512  b=-pp1*pp3*pp4/(p1p2*p2p3*p2p4)
513  c=pp1*pp2*pp4/(p1p3*p2p3*p3p4)
514  d=-pp1*pp2*pp3/(p1p4*p2p4*p3p4)
515  RETURN
516 1 a=0.
517  CALL parab(p2,p3,p4,p,b,c,d)
518  RETURN
519 2 a=0.
520  b=1.
521  c=0.
522  d=0.
523  RETURN
524 3 IF (iabs(k12+k23).eq.2) go to 4
525  go to 1
526 4 d=0.
527  CALL parab(p1,p2,p3,p,a,b,c)
528  RETURN
529  END SUBROUTINE cubic
530  SUBROUTINE parab(p1,p2,p3,p,a,b,c)
531  USE precision
532  IMPLICIT NONE
533  REAL(rprec) :: p1, p2, p3, p, a, b, c, zero, yp, y, xp, x,
534  . p1p2, p2p3, p1p3, p1p4, p2p4, p3p4, pp1, pp2, pp3, pp4
535  INTEGER :: k12, k23, k34
536  data zero/1.e-6/
537 cray lcm (p1),(p2),(p3)
538  p1p2=p1-p2
539  p2p3=p2-p3
540  p1p3=p1-p3
541  pp1=p-p1
542  pp2=p-p2
543  pp3=p-p3
544  k12=sign(1.5_dbl,p1p2)
545  IF(abs(p1p2).lt.zero) go to 1
546  k23=sign(1.5_dbl,p2p3)
547  IF(abs(p2p3).lt.zero) go to 2
548  IF(abs(p1p3).lt.zero) go to 1
549  IF (iabs(k12+k23).ne.2) go to 1
550  a=pp2*pp3/(p1p2*p1p3)
551  b=-pp1*pp3/(p1p2*p2p3)
552  c=pp1*pp2/(p1p3*p2p3)
553  RETURN
554 1 a=0.
555  b=pp3/p2p3
556  c=-pp2/p2p3
557  RETURN
558 2 a=pp3/p1p3
559  b=0.
560  c=-pp1/p1p3
561  RETURN
562  END SUBROUTINE parab
563  SUBROUTINE fit(k,x1,x2,x3,x4,y1,y2,y3,y4,x,y,yp)
564 c --------------------------------------------
565 c this routine is taken from an nyu PROGRAM
566 c set k=1 to find y,yp and k=2 to find x,yp
567 c --------------------------------------------
568  USE precision
569  IMPLICIT NONE
570  REAL(rprec) :: x1, x2, x3, x4, y1, y2, y3, y4, x, y, yp, f,
571  . c1, c2, c3, c4, d1, d2, d3, d4, d12, d13, d14, d23, d24, d34,
572  . crit, xa, xb, ya, yb, dydx
573  INTEGER :: k, iturn, i
574  iturn=0
575  c1=y1/((x1-x2)*(x1-x3)*(x1-x4))
576  c2=y2/((x2-x1)*(x2-x3)*(x2-x4))
577  c3=y3/((x3-x1)*(x3-x2)*(x3-x4))
578  c4=y4/((x4-x1)*(x4-x2)*(x4-x3))
579  IF(k.eq.2) go to 2
580  1 d1=x-x1
581  d2=x-x2
582  d3=x-x3
583  d4=x-x4
584  d12=d1*d2
585  d13=d1*d3
586  d14=d1*d4
587  d23=d2*d3
588  d24=d2*d4
589  d34=d3*d4
590  f=(c1*d23+c2*d13+c3*d12)*d4+c4*d12*d3
591  yp=c1*(d23+d24+d34)+c2*(d13+d14+d34)
592  + +c3*(d12+d14+d24)+c4*(d12+d13+d23)
593  IF(k.eq.2) go to 3
594  y=f
595  RETURN
596  2 IF(y.ge.min(y1,y2,y3,y4).and.y.le.max(y1,y2,y3,y4)) GOTO 4
597  21 CONTINUE
598  WRITE(6,11)x1,x2,x3,x4,y1,y2,y3,y4,y
599  x=0.
600  WRITE(6,500)
601  500 FORMAT("error in sub fit at label 21")
602  WRITE(6,*)'error in furpl--fit1'
603  stop
604  RETURN
605  11 FORMAT(" fit",9e14.5)
606  4 crit=(abs(y1)+abs(y2)+abs(y3)+abs(y4))*1.e-05
607  i=0
608  xa=x1
609  xb=x2
610  ya=y1
611  yb=y2
612  IF((x-xa)*(x-xb).lt.0.) go to 10
613  xa=x2
614  xb=x3
615  ya=y2
616  yb=y3
617  IF((x-xa)*(x-xb).lt.0.) go to 10
618  xa=x3
619  xb=x4
620  ya=y3
621  yb=y4
622  IF((x-xa)*(x-xb).lt.0.) go to 10
623  12 xa=x2
624  ya=y2
625  xb=x3
626  yb=y3
627  x=(xa+xb)/2.
628  IF((y-ya)*(y-yb).lt.0.) go to 1
629  xa=x1
630  ya=y1
631  xb=x2
632  yb=y2
633  x=(xa+xb)/2.
634  IF((y-ya)*(y-yb).lt.0.) go to 1
635  xa=x3
636  ya=y3
637  xb=x4
638  yb=y4
639  x=(xa+xb)/2.
640  IF((y-ya)*(y-yb).lt.0.) go to 1
641  IF(y.ne.y1) go to 13
642  x=x1
643  go to 1
644  13 IF(y.ne.y2) go to 14
645  x=x2
646  go to 1
647  14 IF(y.ne.y3) go to 15
648  x=x3
649  go to 1
650  15 IF(y.ne.y4) go to 16
651  x=x4
652  go to 1
653  16 CONTINUE
654  WRITE(6,11)x,y
655  go to 21
656 3 IF(abs(f-y).lt.crit) iturn=1
657  IF(i.eq.1) go to 7
658  dydx=(yb-ya)/(xb-xa)
659  IF(abs(yp-dydx).lt..2*abs(yp)) go to 7
660  IF((f-y)*(ya-y).lt.0.) go to 5
661  xa=x
662  ya=f
663  go to 6
664  5 xb=x
665  yb=f
666  6 x=(xa+xb)/2.
667  i=1
668  go to 1
669  7 IF((f-y)*(ya-y).lt.0.) go to 8
670  xa=x
671  ya=f
672  go to 9
673  8 xb=x
674  yb=f
675  9 dydx=(yb-ya)/(xb-xa)
676  IF(abs(yp-dydx).lt..2*abs(yp))dydx=yp
677  x=x-(f-y)/dydx
678  IF(iturn.eq.1) RETURN
679  i=0
680  go to 1
681  10 IF((y-ya)*(y-yb).lt.0.) go to 1
682  go to 12
683  END SUBROUTINE fit
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17