ETS  \$Id: Doxyfile 2162 2020-02-26 14:16:09Z g2dpc $
 All Classes Files Functions Variables Pages
scu.f
Go to the documentation of this file.
1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 
3  SUBROUTINE lsq_sur6(r,z,u,n,c,rm,zm,um,dp)
4 c
5 c... list square surf. fitting
6 c c(1)+c(2)*(r-r1)+c(3)*(z-z1)+c(4)*(r-r1)**2+c(5)*(z-z1)**2+c(6)*(r-r1)*(z-z1)
7 c
8  IMPLICIT REAL*8(a-h,o-z)
9 c
10  dimension r(*),z(*),u(*),dp(*)
11  dimension a(6,6),b(6),c(6)
12  dimension ip(6)
13 !!!test DIMENSION u_dum(*)
14 
15  sqrt(xx)=dsqrt(xx)
16 
17 !!!test
18 !
19 ! do i=1,n
20 ! u(i)=
21 ! & 1.d0
22 ! & +2.d0*(r(i)-r(1))
23 ! & +3.d0*(z(i)-z(1))
24 ! & +4.d0*(r(i)-r(1))**2
25 ! & +5.d0*(z(i)-z(1))**2
26 ! & +6.d0*(r(i)-r(1))*(z(i)-z(1))
27 ! enddo
28 !
29 !!!test
30 
31  s_r= 0.d0
32  s_z= 0.d0
33  s_r2=0.d0
34  s_z2=0.d0
35  s_rz=0.d0
36 
37 
38  s_r3= 0.d0
39  s_rz2=0.d0
40  s_r2z=0.d0
41 
42  s_z3=0.d0
43 
44  s_r4= 0.d0
45  s_r2z2=0.d0
46  s_r3z= 0.d0
47 
48  s_z4= 0.d0
49  s_rz3=0.d0
50 
51  s_u= 0.d0
52  s_ur=0.d0
53  s_uz=0.d0
54  s_ur2= 0.d0
55  s_uz2=0.d0
56  s_urz=0.d0
57 
58 
59  do k=1,n
60 
61  s_r= s_r + (r(k)-r(1))
62  s_z= s_z + (z(k)-z(1))
63  s_r2= s_r2 + (r(k)-r(1))**2
64  s_z2= s_z2 + (z(k)-z(1))**2
65  s_rz= s_rz + (r(k)-r(1))*(z(k)-z(1))
66  s_u= s_u +u(k)
67 
68 
69  s_r3= s_r3 + (r(k)-r(1))**3
70  s_rz2= s_rz2 + (r(k)-r(1))*(z(k)-z(1))**2
71  s_r2z= s_r2z + (r(k)-r(1))**2*(z(k)-z(1))
72  s_ur= s_ur + u(k)*(r(k)-r(1))
73 
74  s_z3= s_z3 + (z(k)-z(1))**3
75  s_uz= s_uz + u(k)*(z(k)-z(1))
76 
77  s_r4= s_r4 + (r(k)-r(1))**4
78  s_r2z2= s_r2z2 + (r(k)-r(1))**2*(z(k)-z(1))**2
79  s_r3z= s_r3z + (r(k)-r(1))**3*(z(k)-z(1))
80  s_ur2= s_ur2 + u(k)*(r(k)-r(1))**2
81 
82  s_z4= s_z4 + (z(k)-z(1))**4
83  s_rz3= s_rz3 + (r(k)-r(1))*(z(k)-z(1))**3
84  s_uz2= s_uz2 + u(k)*(z(k)-z(1))**2
85 
86  s_urz= s_urz + u(k)*(r(k)-r(1))*(z(k)-z(1))
87 
88  enddo
89 
90 c the creation the matrix a and right hand b:
91 
92  a(1,1) = dfloat(n)
93  a(1,2) = s_r
94  a(1,3) = s_z
95  a(1,4) = s_r2
96  a(1,5) = s_z2
97  a(1,6) = s_rz
98 
99  b(1) = s_u
100 
101  a(2,1) = s_r
102  a(2,2) = s_r2
103  a(2,3) = s_rz
104  a(2,4) = s_r3
105  a(2,5) = s_rz2
106  a(2,6) = s_r2z
107 
108  b(2) = s_ur
109 
110  a(3,1) = s_z
111  a(3,2) = s_rz
112  a(3,3) = s_z2
113  a(3,4) = s_r2z
114  a(3,5) = s_z3
115  a(3,6) = s_rz2
116 
117  b(3) = s_uz
118 
119  a(4,1) = s_r2
120  a(4,2) = s_r3
121  a(4,3) = s_r2z
122  a(4,4) = s_r4
123  a(4,5) = s_r2z2
124  a(4,6) = s_r3z
125 
126  b(4) = s_ur2
127 
128  a(5,1) = s_z2
129  a(5,2) = s_rz2
130  a(5,3) = s_z3
131  a(5,4) = s_r2z2
132  a(5,5) = s_z4
133  a(5,6) = s_rz3
134 
135  b(5) = s_uz2
136 
137  a(6,1) = s_rz
138  a(6,2) = s_r2z
139  a(6,3) = s_rz2
140  a(6,4) = s_r3z
141  a(6,5) = s_rz3
142  a(6,6) = s_r2z2
143 
144  b(6) = s_urz
145 
146 ! do i=1,6
147 ! amx=0.d0
148 !
149 ! do j=1,6
150 ! if( dabs(a(i,j)) .gt. amx ) then
151 ! amx=dabs(a(i,j))
152 ! endif
153 ! enddo
154 !
155 ! do j=1,6
156 ! a(i,j)=a(i,j)/amx
157 ! enddo
158 !
159 ! b(i)=b(i)/amx
160 !
161 ! enddo
162 
163 c the solution the matrix equation ac=b.
164 
165  CALL ge(6,6,a,b,c,ip)
166 
167 c magnetic axis
168 
169  det=4.d0*c(4)*c(5)-c(6)**2
170 
171  det_r=-2.d0*c(2)*c(5)+c(6)*c(3)
172  det_z=-2.d0*c(3)*c(4)+c(6)*c(2)
173 
174  rm=det_r/det
175  zm=det_z/det
176 
177  um=c(1)+c(2)*rm+c(3)*zm+c(4)*rm**2+c(5)*zm**2+c(6)*rm*zm
178 
179  rm=rm+r(1)
180  zm=zm+z(1)
181 
182  dudr=c(2)
183  dudz=c(3)
184  dudr2=2.d0*c(4)
185  dudz2=2.d0*c(5)
186  dudrdz=c(6)
187 
188  dp(1)=dudr
189  dp(2)=dudz
190  dp(3)=dudr2
191  dp(4)=dudrdz
192  dp(5)=dudz2
193 
194 
195  RETURN
196  END
197 
198 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
199 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200 
201  SUBROUTINE deriv5(X,Y,F,M,N,U)
202 c
203 c... definiton of the first and second derivations.
204 c
205 c
206  IMPLICIT REAL*8(a-h,o-z)
207 c
208  dimension x(*),y(*),f(*),u(*)
209  dimension a(5,5),b(5),s(5)
210 c dimension aa(5,5),w1(5),w2(5)
211  dimension ip(5)
212 c
213  sqrt(r)=dsqrt(r)
214 c
215 c
216  sdx2 = 0.d0
217  sdy2 = 0.d0
218  sdxdy = 0.d0
219  sdx3 = 0.d0
220  sdx2dy= 0.d0
221  sdxdy2= 0.d0
222  sdy3 = 0.d0
223  sdx4 = 0.d0
224  sdx3dy= 0.d0
225  sdx2y2= 0.d0
226  sdxdy3= 0.d0
227  sdy4 = 0.d0
228 c
229  daver2= 0.d0
230 c
231  fdx = 0.d0
232  fdy = 0.d0
233  fdx2 = 0.d0
234  fdxdy = 0.d0
235  fdy2 = 0.d0
236 c
237 c the main loop
238 c
239  DO 1 i=2,m
240  dx=x(i)-x(1)
241  dy=y(i)-y(1)
242  df=f(i)-f(1)
243  dx2=dx*dx
244  dy2=dy*dy
245 c
246  sdx2 =sdx2 + dx2
247  sdy2 =sdy2 + dy2
248  sdxdy =sdxdy + dx *dy
249  sdx3 =sdx3 + dx2 *dx
250  sdx2dy=sdx2dy + dx2 *dy
251  sdxdy2=sdxdy2 + dx *dy2
252  sdy3 =sdy3 + dy2 *dy
253  sdx4 =sdx4 + dx2 *dx2
254  sdx3dy=sdx3dy + dx2 *dx *dy
255  sdx2y2=sdx2y2 + dx2 *dy2
256  sdxdy3=sdxdy3 + dx *dy2 *dy
257  sdy4 =sdy4 + dy2 *dy2
258 c
259  daver2=daver2 + dx2 + dy2
260 c
261  fdx =fdx + df *dx
262  fdy =fdy + df *dy
263  fdx2 =fdx2 + df *dx2
264  fdxdy =fdxdy + df *dx *dy
265  fdy2 =fdy2 + df *dy2
266 c
267  1 CONTINUE
268 c
269 c the creation the matrix a:
270 c
271  a(1,1) = sdx2
272  a(1,2) = sdxdy
273  a(1,3) = sdx3 * 0.5d0
274  a(1,4) = sdx2dy
275  a(1,5) = sdxdy2 * 0.5d0
276  a(2,1) = sdxdy
277  a(2,2) = sdy2
278  a(2,3) = sdx2dy * 0.5d0
279  a(2,4) = sdxdy2
280  a(2,5) = sdy3 * 0.5d0
281  a(3,1) = sdx3
282  a(3,2) = sdx2dy
283  a(3,3) = sdx4 * 0.5d0
284  a(3,4) = sdx3dy
285  a(3,5) = sdx2y2 * 0.5d0
286  a(4,1) = sdx2dy
287  a(4,2) = sdxdy2
288  a(4,3) = sdx3dy * 0.5d0
289  a(4,4) = sdx2y2
290  a(4,5) = sdxdy3 * 0.5d0
291  a(5,1) = sdxdy2
292  a(5,2) = sdy3
293  a(5,3) = sdx2y2 * 0.5d0
294  a(5,4) = sdxdy3
295  a(5,5) = sdy4 * 0.5d0
296 c
297 c the creation the right hand b:
298 c
299  b(1) = fdx
300  b(2) = fdy
301  b(3) = fdx2
302  b(4) = fdxdy
303  b(5) = fdy2
304 c
305 c the weights s(i):
306 c
307  daver=sqrt(daver2)
308 c
309  s(1) = 1.d0/sqrt(sdx2)/daver
310  s(2) = 1.d0/sqrt(sdy2)/daver
311  s(3) = 1.d0/sdx2/daver
312  s(4) = 1.d0/sqrt(sdx2*sdy2)/daver
313  s(5) = 1.d0/sdy2/daver
314 c
315  DO 2 i=1,n
316  DO 3 j=1,n
317  a(i,j)=s(i)*a(i,j)
318  3 CONTINUE
319  b(i)=s(i)*b(i)
320  2 CONTINUE
321 c
322 c the solution the matrix equation au=b.
323 c
324  ifail=0
325 c
326 ccccc CALL f04atf(a,5,b,n,u,aa,5,w1,w2,ifail)
327  CALL ge(5,5,a,b,u,ip)
328 c
329 ccccc IF( ifail .NE. 0 ) WRITE (6,*)' IFAIL= ',ifail
330 c
331  RETURN
332  END
333 
334 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
335 c
336 c... cubic fit
337 c
338  SUBROUTINE deriv9(X,Y,F,M,N,U)
339 c
340  IMPLICIT REAL*8(a-h,o-z)
341 c
342  dimension x(1:m),y(1:m),f(1:m),u(1:n)
343  dimension a(1:9,1:9),b(1:9)
344  dimension ip(1:9)
345 ! DIMENSION AA(1:100,1:9),BB(1:100)
346  REAL*8, ALLOCATABLE :: aa(:,:),bb(:)
347 c
348 c... dimension check
349 ! IF(M.GT.100) THEN
350 ! WRITE(6,*) ' DERIV9: M EXCEEDS 100'
351 ! STOP
352 ! ENDIF
353  IF (.NOT. ALLOCATED(aa)) ALLOCATE(aa(m,9))
354  IF (.NOT. ALLOCATED(bb)) ALLOCATE(bb(m))
355 c... the main loop
356 c
357  DO i=1,m-1
358  dx=x(i+1)-x(1)
359  dy=y(i+1)-y(1)
360  df=f(i+1)-f(1)
361  dxx=.5*dx*dx
362  dxy= dx*dy
363  dyy=.5*dy*dy
364 c
365  aa(i,1)=dx
366  aa(i,2)=dy
367  aa(i,3)=dxx
368  aa(i,4)=dxy
369  aa(i,5)=dyy
370  !cubic
371  dxxx= dx*dx*dx
372  dxxy= dx*dx*dy
373  dxyy= dx*dy*dy
374  dyyy= dy*dy*dy
375  aa(i,6)=dxxx
376  aa(i,7)=dxxy
377  aa(i,8)=dxyy
378  aa(i,9)=dyyy
379 c
380  bb(i)=df
381 c
382  ENDDO
383 c
384 c... the creation the matrix a:
385 c
386  DO k=1,9
387  DO l=1,9
388  akl=0.
389  DO i=1,m-1
390  akl=akl+aa(i,k)*aa(i,l)
391  ENDDO
392  a(k,l)=akl
393  ENDDO
394  ENDDO
395 c
396 c... the creation the right hand b:
397 c
398  DO k=1,9
399  bk=0.
400  DO i=1,m-1
401  bk=bk+aa(i,k)*bb(i)
402  ENDDO
403  b(k)=bk
404  ENDDO
405 c
406  DEALLOCATE(aa,bb)
407 c
408 c... the solution the matrix equation au=b.
409 c
410  ifail=0
411 c
412  CALL ge(n,9,a,b,u,ip)
413 c
414  IF( ifail .NE. 0 ) WRITE (6,*)' DERIV9:IFAIL= ',ifail
415 c
416  RETURN
417  END
418 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
419 
420  SUBROUTINE ge(N,NZ,A,X,Y,IP)
421 c
422 c... gauss elimination
423 c
424  IMPLICIT REAL*8(a-h,o-z)
425 c
426  dimension a(1:nz,1:nz),x(1:n),y(1:n),ip(1:n)
427 c
428  abs(xarg)=dabs(xarg)
429 c
430  DO 90 i=1,n
431  ip(i)=i
432  90 CONTINUE
433 c
434  DO 100 i=1,n-1
435 c... max row element search
436  rm=abs(a(i,i))
437  jm=i
438  DO 101 j=i,n
439  aba=abs(a(i,j))
440  IF(aba.GT.rm) THEN
441  rm=aba
442  jm=j
443  ENDIF
444  101 CONTINUE
445 c... permutations
446  ipe=ip(i)
447  ip(i)=ip(jm)
448  ip(jm)=ipe
449  DO 102 k=1,n
450  ape=a(k,i)
451  a(k,i)=a(k,jm)
452  a(k,jm)=ape
453  102 CONTINUE
454  ad=1.d0/a(i,i)
455  DO 100 k=i+1,n
456  am=a(k,i)*ad
457  x(k)=x(k)-am*x(i)
458  DO 100 j=i,n
459  a(k,j)=a(k,j)-am*a(i,j)
460  100 CONTINUE
461 c
462  y(n)=x(n)/a(n,n)
463  DO 200 i=n-1,1,-1
464  ad=1./a(i,i)
465  y(i)=x(i)
466  DO 201 j=n,i+1,-1
467  y(i)=y(i)-a(i,j)*y(j)
468  201 CONTINUE
469  y(i)=y(i)*ad
470  200 CONTINUE
471 c
472 c... back permutation
473  DO 300 i=1,n
474  x(ip(i))=y(i)
475  300 CONTINUE
476  DO 400 i=1,n
477  y(i)=x(i)
478  400 CONTINUE
479 c
480  RETURN
481  END
482 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483 c
484  REAL*8 FUNCTION greeni(R,Z,RP,ZP)
485 c green'S FUNCTION FOR TOROIDAL CURRENT LOOP IN INFINITY AREA.
486 C
487  IMPLICIT REAL*8(A-H,O-Z)
488 C
489  REAL*8 MMDELK,MMDELE
490 C
491  SQRT(X)=DSQRT(X)
492 C
493  T=SQRT(4.d0*R*RP/((R+RP)**2+(Z-ZP)**2))
494 C
495 C CALCULATION OF COMPLETE ELLIPTIC INTEGRALS OF 1TH AND 2TH KIND.
496 C
497  TT=(1.d0-T**2)
498 C
499 C THE FUNCTIONS MMDELK, MMDELE FROM IMSL LIBRARY.
500 C
501 C ELCK=MMDELK(2,T,IFAILK)
502 C ELCE=MMDELE(2,T,IFAILE)
503 C
504 C THE FUNCTIONS S21BBF, S21BCF FROM NAG LIBRARY.
505 C
506  ELCK= S21BBF(0.D0,TT,1.D0,IFAILK)
507  ELCE=ELCK-T*T/3.D0*S21BCF(0.D0,TT,1.D0,IFAILE)
508 C
509 .NE..OR..NE. IF(IFAILK0IFAILE0) WRITE(*,*)' attention! GREEN F.:',
510  , ' IFAILK, IFAILE= ',ifailk,ifaile
511 c
512  greeni = ( (1.d0-t*t*0.5d0)*elck-elce )*( sqrt(r*rp)/t )
513 c
514  RETURN
515  END
516 c
517 c********************************************************************
518 
519  subroutine grgren(R0,Z0,R,Z, dGdr,dGdz)
520 
521 c green's function derivatives at the point (R,Z)
522 c the source position is (R0,Z0)
523 c to obtain the real values of derivatives dGdr,dGdz must be devided by pi
524 
525  IMPLICIT REAL*8(A-H,O-Z)
526  SQRT(X)=DSQRT(X)
527 
528  t=SQRT( 4.D0*R*R0/( (R+R0)**2+(Z-Z0)**2 ) )
529 
530  tt=(1.D0-t**2)
531 
532 C THE FUNCTIONS S21BBF, S21BCF FROM NAG LIBRARY.
533 
534  fK= S21BBF(0.D0,tt,1.D0,IFAILK)
535  fE=fK-t*t/3.D0*S21BCF(0.D0,tt,1.D0,IFAILE)
536 
537 .NE..OR..NE. IF(IFAILK0IFAILE0) WRITE(*,*)' attention! GREEN F.:',
538  * ' IFAILK, IFAILE= ',ifailk,ifaile
539 
540  q = ( (1.d0-t*t*0.5d0)*fk-fe )/t
541 
542  dtdr = 0.5d0*( t/r-t*t*t*(r+r0)/(2.d0*r*r0) )
543 
544  dtdz =-t*t*t*(z-z0)/(4.d0*r*r0)
545 
546  dqdt = -q/t + 0.5d0*( fe/tt-fk)
547 
548  dgdr = 0.5d0*sqrt(r0/r)*q + sqrt(r0*r)*dqdt*dtdr
549 
550  dgdz = sqrt(r0*r)*dqdt*dtdz
551 
552  RETURN
553  END
554 
555 c********************************************************************
556 
557  SUBROUTINE bint(X,Y,R0,Z0,r1,z1,F,I)
558  IMPLICIT REAL*8(a-h,o-z)
559 c
560 c===============functions===============functions================c
561 c
562  sqrt(x)=dsqrt(x)
563  alog(x)=dlog(x)
564 c
565 c======================special functions=========================c
566 c
567  d(x1,y1,x2,y2)= sqrt( (x1-x2)**2 + (y1-y2)**2 )
568 c
569  tint(t)=(t-pscal)*(0.5d0*dlog(pvec**2+(t-pscal)**2+eps**2)-1.d0)+
570  + pvec*datan((t-pscal)/(pvec+eps))
571 c
572 c===============functions===============functions================c
573 c
574  epsh=1.d-9
575  f=0.d0
576 
577  w=0.d0
578  h=d(r0,z0,r1,z1)
579  eps=epsh*h
580 c
581 c integral along the edge of region
582 c
583  rm=0.5d0*(r0+r1)
584  zm=0.5d0*(z0+z1)
585  dtm=d(x,y,rm,zm)
586  dfm=d(-x,y,rm,zm)
587 c
588  IF( i .NE. 0 ) THEN
589 c
590  pscal=( (-x-r0)*(r1-r0) + (y-z0)*(z1-z0) )/h
591  pvec= ( (-x-r0)*(z1-z0) - (y-z0)*(r1-r0) )/h
592 c
593  w = tint(h)-tint(0.d0)
594 c
595  pscal=( (x-r0)*(r1-r0) + (y-z0)*(z1-z0) )/h
596  pvec= ( (x-r0)*(z1-z0) - (y-z0)*(r1-r0) )/h
597 c
598  w = w - (tint(h)-tint(0.d0))
599 c
600  IF( dtm .LT. 0.25d0*h ) THEN
601  dfj =d(-x,y,r0,z0)
602  dfj1=d(-x,y,r1,z1)
603  dtj =d( x,y,r0,z0)
604  dtj1=d( x,y,r1,z1)
605 c
606 c+++ w = w - h * alog( dfj*dfj1/(d(x,y,r0,z0)*d(x,y,r1,z1)) )
607 c
608  w = w*dfm-0.5d0*h*(alog(dfj/dtj)*dfj+alog(dfj1/dtj1)*dfj1)
609 c
610  ELSE
611 c
612 c+++ w = w - h * alog( dfm / dtm )
613 c
614  w = ( w - h * alog( dfm / dtm ) ) * dfm
615 c
616  ENDIF
617  ENDIF
618 c
619  IF( dtm .LT. 0.25d0*h ) THEN
620 c
621  e = 0.5d0 * ( greeni(x,y,r0,z0) + greeni(x,y,r1,z1) )
622 c
623  ELSE
624 c
625  e = greeni(x,y,rm,zm)
626 c
627  ENDIF
628 c
629  f = e*h + 0.25d0*w
630 c
631  f = f /3.14159265359d0
632 c
633  RETURN
634  END
635 
636 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
637  SUBROUTINE bint_d(X,Y,R0,Z0,r1,z1,Fr,Fz,i)
638 
639  IMPLICIT REAL*8(a-h,o-z)
640 
641  sqrt(xxx)=dsqrt(xxx)
642 
643  d(x1,y1,x2,y2)= dsqrt( (x1-x2)**2 + (y1-y2)**2 )
644 
645  rm=0.5d0*(r0+r1)
646  zm=0.5d0*(z0+z1)
647 
648  call grgren(rm,zm,x,y, dgdr,dgdz)
649 
650  dell=d(r0,z0,r1,z1)
651 
652  fr=dgdr*dell/3.14159265359d0
653  fz=dgdz*dell/3.14159265359d0
654 
655  RETURN
656  END
657 c********************************************************************
658  subroutine greng(R0,Z0,R,Z,fgreen,dGdr,dGdz)
659 
660 c green's function derivatives at the point (R,Z)
661 c the source position is (R0,Z0)
662 c to obtain the real values of derivatives dGdr,dGdz must be devided by pi
663 
664  IMPLICIT REAL*8(A-H,O-Z)
665  SQRT(X)=DSQRT(X)
666 
667  t=SQRT( 4.D0*R*R0/( (R+R0)**2+(Z-Z0)**2 ) )
668 
669  tt=(1.D0-t**2)
670 
671 C THE FUNCTIONS S21BBF, S21BCF FROM NAG LIBRARY.
672 
673  fK= S21BBF(0.D0,tt,1.D0,IFAILK)
674  fE=fK-t*t/3.D0*S21BCF(0.D0,tt,1.D0,IFAILE)
675 
676 .NE..OR..NE. IF(IFAILK0IFAILE0) WRITE(*,*)' attention! GREEN F.:',
677  * ' IFAILK, IFAILE= ',ifailk,ifaile
678 
679  q = ( (1.d0-t*t*0.5d0)*fk-fe )/t
680 
681  fgreen= q*sqrt(r*r0)
682 
683  dtdr = 0.5d0*( t/r-t*t*t*(r+r0)/(2.d0*r*r0) )
684 
685  dtdz =-t*t*t*(z-z0)/(4.d0*r*r0)
686 
687  dqdt = -q/t + 0.5d0*( fe/tt-fk)
688 
689  dgdr = 0.5d0*sqrt(r0/r)*q + sqrt(r0*r)*dqdt*dtdr
690 
691  dgdz = sqrt(r0*r)*dqdt*dtdz
692 
693  RETURN
694  END
695 c********************************************************************
696  subroutine d2gren(R0,Z0,R,Z,d2Gdrz,d2Gdzz,d2Gdrr)
697 
698 c green's function second derivatives at the point (R,Z)
699 c the source position is (R0,Z0)
700 c to obtain the real values of derivatives dGdr,dGdz must be devided by pi
701 
702  IMPLICIT REAL*8(A-H,O-Z)
703  SQRT(X)=DSQRT(X)
704 
705  t=SQRT( 4.D0*R*R0/( (R+R0)**2+(Z-Z0)**2 ) )
706 
707  tt=(1.D0-t**2)
708 
709 C THE FUNCTIONS S21BBF, S21BCF FROM NAG LIBRARY.
710 
711  fK= S21BBF(0.D0,tt,1.D0,IFAILK)
712  fE=fK-t*t/3.D0*S21BCF(0.D0,tt,1.D0,IFAILE)
713 
714 .NE..OR..NE. IF(IFAILK0IFAILE0) WRITE(*,*)' attention! GREEN F.:',
715  * ' IFAILK, IFAILE= ',ifailk,ifaile
716 
717  q = ( (1.d0-t*t*0.5d0)*fk-fe )/t
718 
719  fgreen= q*sqrt(r*r0)
720 
721  dtdr = 0.5d0*( t/r-t*t*t*(r+r0)/(2.d0*r*r0) )
722 
723  dtdz =-t*t*t*(z-z0)/(4.d0*r*r0)
724 
725  dqdt = -q/t + 0.5d0*( fe/tt-fk)
726 
727  d2tdrz =
728  & -(1/r/2-3.d0/4.d0*t**2*(r+r0)/r/r0)*t**3*(z-z0)/r/r0/4.d0
729 
730 
731  d2tdzz =
732  & 3.d0/16.d0*t**5*(z-z0)**2/r**2/r0**2-t**3/r/r0/4.d0
733 
734 
735  d2tdrr=
736  &(1/r/2-3.d0/4.d0*t**2*(r+r0)/r/r0)*(t/r/2-t**3*(r+r0)/r/r0/4)
737  #-t/r**2/2-t**3/r/r0/4+t**3*(r+r0)/r**2/r0/4
738 
739 
740  d2qdt2 =
741  & q/t**2+fe/(1-t**2)**2*t+1/(1-t**2)*(fe-fk)/t/2-(fe/(1-t**2)-fk)
742  #/t/2.D0
743 
744 
745  dgdr = 0.5d0*sqrt(r0/r)*q + sqrt(r0*r)*dqdt*dtdr
746  d2gdzz=sqrt(r*r0)*(d2qdt2*dtdz**2+dqdt*d2tdzz)
747  d2gdrz=sqrt(r*r0)*(dqdt*dtdz/r/2.d0+d2qdt2*dtdr*dtdz+dqdt*d2tdrz)
748  d2gdrr=sqrt(r*r0)*(dqdt*dtdr/r/2.d0+d2qdt2*dtdr*dtdr+dqdt*d2tdrr-
749  & q/r**2/2.d0) + dgdr/r/2.d0
750 
751  !dGdr = 0.5d0*sqrt(R0/R)*Q + sqrt(R0*R)*dQdt*dtdr
752  !dGdz = sqrt(R0*R)*dQdt*dtdz
753 
754 
755  RETURN
756  END
757 
758 c********************************************************************
759 
760  REAL*8 FUNCTION cpr(H,NP,S)
761 
762  IMPLICIT REAL*8(a-h,o-z)
763 
764  f(x)=h*(1.d0-x**np)/(1.d0-x)-s
765 
766  i=0
767  x=1.d0
768  eps=1.d-8
769 
770  a=0.0001d0
771  b=20.d0
772 
773  3 x=(a+b)/2.d0
774  y0=f(x)
775  y1=f(a)
776  y2=f(b)
777  i=i+1
778 
779  if(i.gt.1000) then
780 
781  write(6,*) ' ***print from routine CPR:'
782  write(6,*) ' can`t generate grid out of plasma domain'
783  write(6,*) ' program is terminated'
784  stop
785 
786  endif
787 
788  IF(dabs(y0).LT.eps) go to 1
789  IF((y0*y2).GT.0.d0) go to 2
790 
791  a=x
792  go to 3
793  2 b=x
794  go to 3
795  1 CONTINUE
796  cpr=x
797 
798  RETURN
799  END
800 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
801 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
802 
803  REAL*8 FUNCTION funsq(R1,R2,R3,R4,Z1,Z2,Z3,Z4)
804 
805  include 'double.inc'
806 
807  r13=r3-r1
808  r24=r4-r2
809  z13=z3-z1
810  z24=z4-z2
811  funsq=(r13*z24-r24*z13)*0.5d0
812 
813  !IF(ZS.LT.0.) WRITE(6,*) '***ERROR:SQUARE<0',zs
814 
815  !FUNSQ=ZS*0.5D0
816 
817  RETURN
818  END
819 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
820 
821 
822  SUBROUTINE prog1d(M,U, A,B,C,F, ALF,BET)
823 !**********************************************************************
824 ! PROGONKA - USUAL
825 !
826 ! A(J)*U(J+1)-C(J)*U(J)+B(J)*U(J-1)=-F(J) , J=1,M
827 !
828 ! B(1)=0 , A(M)=0
829 !
830 !**********************************************************************
831 
832  IMPLICIT REAL*8(a-h,o-z)
833 
834  dimension u(*), a(*),b(*),c(*),f(*), alf(*),bet(*)
835 
836  alf(1)=a(1)/c(1)
837  bet(1)=f(1)/c(1)
838 
839  DO 20 k=2,m
840 
841  dk=c(k)-alf(k-1)*b(k)
842  alf(k)=a(k)/dk
843 20 bet(k)=(f(k)+bet(k-1)*b(k))/dk
844 
845  u(m) = bet(m)
846 
847  DO 30 j=2,m
848 
849  k=m+1-j
850 30 u(k)=alf(k)*u(k+1)+bet(k)
851 
852  RETURN
853  END
854 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
855 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
856 
857  real*8 function blin_(r0,z0,r1,r2,z1,z2,u1,u2,u3,u4 )
858 
859  include 'double.inc'
860 
861  s1=(r2-r0)*(z2-z0)
862  s3=(r0-r1)*(z0-z1)
863  s2=(r0-r1)*(z2-z0)
864  s4=(r2-r0)*(z0-z1)
865 
866  u0=(s1*u1+s2*u2+s3*u3+s4*u4)/(s1+s2+s3+s4)
867 
868  blin_=u0
869 
870  return
871  end
872 
873 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
874 
875 
876 
877 
878 
879 
funsq
REAL *8 function funsq(R1, R2, R3, R4, Z1, Z2, Z3, Z4)
Definition: scu.f:803
bint_d
subroutine bint_d(X, Y, R0, Z0, r1, z1, Fr, Fz, i)
Definition: scu.f:637
greeni
REAL *8 function greeni(R, Z, RP, ZP)
Definition: scu.f:484
prog1d
subroutine prog1d(M, U, A, B, C, F, ALF, BET)
Definition: scu.f:822
deriv9
subroutine deriv9(X, Y, F, M, N, U)
Definition: scu.f:338
bint
subroutine bint(X, Y, R0, Z0, r1, z1, F, I)
Definition: scu.f:557
deriv5
subroutine deriv5(X, Y, F, M, N, U)
Definition: scu.f:201
analytical_functions::ad
REAL(R8) function ad(IION, X, T)
Definition: analytical_functions.f90:281
loop
subroutine loop
Definition: Eq2_m.f:231
axis
subroutine axis(n, h, r, f)
Definition: solution2.f90:586
blin_
real *8 function blin_(r0, z0, r1, r2, z1, z2, u1, u2, u3, u4)
Definition: scu.f:857
lsq_sur6
subroutine lsq_sur6(r, z, u, n, c, rm, zm, um, dp)
Definition: scu.f:3
phys_functions::r2
real(r8) function r2(a, x, xr, xs, yr, ys, psi, psir, F_dia)
Definition: phys_functions.f90:121
ge
subroutine ge(N, NZ, A, X, Y, IP)
Definition: scu.f:420
grgren
subroutine grgren(R0, Z0, R, Z, dGdr, dGdz)
Definition: scu.f:519
d2gren
subroutine d2gren(R0, Z0, R, Z, d2Gdrz, d2Gdzz, d2Gdrr)
Definition: scu.f:696
integral
subroutine integral(n, h, r, f, int)
Definition: solution2.f90:527
greng
subroutine greng(R0, Z0, R, Z, fgreen, dGdr, dGdz)
Definition: scu.f:658
cpr
REAL *8 function cpr(H, NP, S)
Definition: scu.f:760
matrix
subroutine matrix
Definition: Matrix.f:2