V3FIT
garc.f
1  SUBROUTINE garc(tp,xp,zp,csx,csz,npmax,arc,npc)
2  USE precision
3  IMPLICIT NONE
4  INTEGER :: npmax,npc
5  REAL(rprec) :: tp(*),xp(*),zp(*),arc(*)
6  REAL(rprec) :: csx(*),csz(*)
7 c REAL(rprec) :: csx(3,npmax),csz(3,npmax)
8  REAL(rprec) :: yg(4),wg(4)
9  REAL(rprec) :: el,summ,a,b
10  INTEGER ngaus,j,is,i
11  REAL(rprec) :: ws1,ws2,ws,delta
12 c ----------------------------------------------------------
13 c gaussian quadrature constants for the INTerval from 0 to 1
14 c ----------------------------------------------------------
15  data yg/.069431844202974,.330009478207572,.669990521792428,
16  1 .930568155797027/,wg/.173927422568727,
17  2 2*.326072577431273,.173927422568727/,ngaus/4/
18  arc(1)=0.
19  el=0.
20  DO j=2,npc
21  is=j-1
22 c --------------------------------------------------------------------
23 c this routine USEs a four poINT gaussian quadrature to compute the
24 c INTegral of SQRT((dx/dt)**2+(dz/dt)**2) with respect to t around the
25 c curve. the quadrature is on poINT number tp.
26 c --------------------------------------------------------------------
27  summ=0.
28  delta=tp(j)-tp(is) ! delta is h
29  DO i=1,ngaus
30  b=yg(i)
31  a=1.-yg(i)
32  ws1=(xp(j)-xp(is))/delta-(3.*a**2-1.)/6.*delta*csx(is)
33  $ +(3.*b**2-1.)/6.*delta*csx(j)
34  ws2=(zp(j)-zp(is))/delta-(3.*a**2-1.)/6.*delta*csz(is)
35  $ +(3.*b**2-1.)/6.*delta*csz(j)
36 c yq=yg(i)*delta !yq is a*h, yg(i) is b,
37 c ws1=(3.*csx(3,is)*yq+2.*csx(2,is))*yq+csx(1,is)
38 c ws2=(3.*csz(3,is)*yq+2.*csz(2,is))*yq+csz(1,is)
39  ws=sqrt(ws1*ws1+ws2*ws2)
40  summ=summ+ws*wg(i)
41  END DO
42  el=el+summ*delta
43  arc(j)=el
44  END DO
45  RETURN
46  END