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