Editing 2DStructure/ToroidalCoordinateIntegrationLimits (section)

==Coming Up to Speed==
I am returning to this problem after letting it lay dormant for more than 14 months.  These notes are being recorded as I try to decipher the fortran code that I have previously developed.  Here is a screenshot from philip.hpc.lsu.edu that will help me archive the dates various fortran routines were last edited.
<div align="center">
[[File:PhilipScreenshot01.png|400px|Screenshot of philip.hpc.lsu.edu]][[File:WarpedSurfaceZ0.png|400px|Added Z0 = 0 seam to warped surface]]
<!--[[File:PlayWithPieces.png|350px|Early playing with ZoneII pieces]]-->
</div>
Let's start by reviewing the pair of routines (''mainI_03.for'',''volumeI_03.for'') which, presumably, integrate over the region labeled, [[#Zone_I:|Zone I]], above.

===Main Routine===

<div id="MainRoutine" style="width: 85%; height: 50em; overflow: auto;">
<table border="1" align="center" width="90%">
<tr>
  <th align="center">Unedited ''mainI_03.for''</th>
  <th align="center">Current Comments &amp; Notes</th>
</tr>
<tr>
  <td align="left">
<pre>
June 12, 2016
ZoneI -- mainI_03.for
(Generates good Vistrails error map.)
========================================

PROGRAM MainI
      real*8 dxx,dzz,dvarpi_t,dr_t,TorusPot,Volume,total,error
      real*8 p3,p5
      real*8 griddx,Xstart,Xend,Ystart,Yend,gridedge
      real*8 x(49),y(49),xfull(50),yfull(50)
      real*8 e(49,49),esum
      real*8 elog
      real   emax,emin
      real   esumshort,elimit
      real   eshort(49,49),xshort(50),yshort(50)
      integer iiZone,ngrid,i,j,k,ncount
  151 format(1x,'Zone I:',5x,'Xstart, Xend, Ystart, Yend = ',&
     &1P4d14.5)
  152 format(1x,' j ',' k ',7X,'x',14X,'y',10X,'error')
  153 format(2I5,1p3d15.6)
  154 format(5x,'Across ',I6,' zones, the average error = ',1p1d15.2)
  601 format(1x,'emax, emin = ',1p2E15.3)
      dvarpi_t = 0.75d0
      dr_t = 0.25d0
      ngrid = 50
      gridedge = 1.5d0
      griddx = gridedge/dfloat(ngrid-1)
      xfull(1)=0.0d0 
      yfull(1)=0.0d0 
      do i=2,ngrid
      xfull(i) = xfull(i-1)+griddx
      yfull(i) = yfull(i-1)+0.5d0*griddx
      enddo
      x(1)=0.5d0*griddx 
      y(1)=0.5d0*griddx 
      do i=2,ngrid-1
      x(i) = x(i-1)+griddx
      y(i) = y(i-1)+0.5d0*griddx
      enddo
      do j=1,ngrid-1
        do k=1,ngrid-1
        e(j,k)=0.0d0
        enddo
      enddo
!******
!
!  Zone I
!
!******
! [iiZone = 0] Full Volume
!
      iizone = 0
      Xstart = 0.0d0
      Xend   = gridedge
      Ystart = dr_t
      Yend   = y(ngrid-1)+0.5d0*griddx
      write(*,151)Xstart,Xend,Ystart,Yend
      write(*,152)
      ncount=0
      esum=0.0d0
      do k=1,ngrid-1
        do j=1,ngrid-1
          if(x(j).gt.Xstart .and. x(j).lt.Xend)then
            if(y(k).gt.Ystart .and. y(k).lt.Yend)then
            dxx = x(j)
            dzz = y(k)
        total = 0.0d0
        call Integrate(dxx,dzz,dvarpi_t,dr_t,iiZone,TorusPot,Volume,p3,p5)
        e(j,k) = (1.0d0-Volume)
        ncount=ncount+1
        esum = esum + DABS(e(j,k))
        write(*,153)j,k,x(j),y(k),e(j,k)
            endif
          endif
        enddo
      enddo
      esum = esum/ncount
      write(*,154)ncount,esum
! Prepare for XMLwriter
      do k=1,ngrid
      xshort(k) = xfull(k)
      yshort(k) = yfull(k)
      enddo
      do k=1,ngrid-1
        do j=1,ngrid-1
        if(e(j,k).eq.0.0d0)e(j,k)=esum
        enddo
      enddo
      do k=1,ngrid-1
        do j=1,ngrid-1
        elog = dlog10(DABS(e(j,k)))
        eshort(j,k) = DABS(elog)
        enddo
      enddo
      emax = 0.0
      emin = 20.0
      esumshort = esum
      elimit = 1.5*esumshort
      do k=1,ngrid-1
        do j=1,ngrid-1
        if(eshort(j,k).gt.emax)emax=eshort(j,k)
        if(eshort(j,k).lt. emin)emin=eshort(j,k)
        if(eshort(j,k).lt.elimit)eshort(j,k)=elimit
        enddo
      enddo
      write(*,601)emax,emin
      do k=1,ngrid-1
        do j=1,ngrid-1
        eshort(j,k) = (eshort(j,k)-elimit)/(emax-elimit)
        enddo
      enddo
      emax = 0.0
      emin = 20.0
      do k=1,ngrid-1
        do j=1,ngrid-1
        if(eshort(j,k).gt.emax)emax=eshort(j,k)
        if(eshort(j,k).lt. emin)emin=eshort(j,k)
        enddo
      enddo
      write(*,601)emax,emin
      do k=1,ngrid-1
        do j=1,ngrid-1
        eshort(j,k) = (eshort(j,k)-emin)/(emax-emin)
        enddo
      enddo
      call XMLwriter01(ngrid,xshort,yshort,eshort)
      stop
END PROGRAM MainI
</pre>
  </td>
  <td align="left" width="50%">
<ul>
<li>dvarpi_t: &nbsp; <math>~\varpi_t = \tfrac{3}{4}</math></li>
<li>dr_t: &nbsp; <math>~r_t = \tfrac{1}{4}</math></li>
<li>(xfill, yfill) marks vertices, 1 --> ngrid = 50</li>
<li>(x, y)  marks cell centers, 1 --> (ngrid-1)</li>
<li>Rectilinear grid with <math>~\Delta x = \Delta y = 1.5/(50 - 1)</math></li>
<li>iiZone = 0</li>
<ul>
<li><font color="red">(Xstart, Xend, Ystart, Yend)</font> = (0, 1.5, r_t, 1.5)</li>
<li>(nested double) Loop through all cell ''centers'' that fall within these just-defined (X, Y) boundaries</li>
<li>For each cell center:</li>
<ul>
<li>call ''Integrate''</li>
<li>calculate the "error", e = (1 - Volume)</li>
<li>update the "counter", ncount = ncount + 1
<li>add absolute value of the "error" to accumulated error, esum
<li>write: <font color="red">j, k, x, y, error</font></li>
</ul>
<li>esum = esum/ncount</li>
<li>write: <font color="red">ncount, esum</font></li>
</ul>
</ul>


----


<ul>
  <li>Arguments of Integrate(dxx,dzz,dvarpi_t,dr_t,iiZone,TorusPot,Volume,p3,p5)</li>
  <ul>
    <li>dxx: &nbsp; cell-centered, x(j)</li>
    <li>dzz: &nbsp; cell-centered, y(k)</li>
    <li><math>~\varpi_t</math></li>
    <li><math>~r_t</math></li>
    <li>iiZone</li>
    <li><font color="darkblue">Return:</font> &nbsp; TorusPot</li>
    <li><font color="darkblue">Return:</font> &nbsp; Volume</li>
    <li>p3</li>
    <li>p5</li>
  </ul>
</ul>
  </td>
</tr>
</table>
</div>
&nbsp;<br />
<font color="red"><b>Attention!</b></font> &nbsp; As a default, only a portion of the code is displayed in the above table; using your mouse, point anywhere inside the table then scroll down/up to view the entire file.
&nbsp;<br />
&nbsp;<br />

===Testing Volume Integration===

<div id="TestingVolumeIntegration" style="width: 85%; height: 50em; overflow: auto;">

<table border="1" align="center" width="90%">
<tr>
  <th align="center">Unedited ''volumeI_03.for''</th>
  <th align="center">Current Comments &amp; Notes</th>
</tr>
<tr>
  <td align="left">
<pre>
      subroutine Integrate(RR,ZZ,v_t,r_t,iiZone,TorusPot,sumVol,pieceIII3,pieceIII5)
!!!!!!!!!!!!
!
!  Validate PATTERN III-B
!
!!!!!!!!!!!!
      real*8 TorusPot,RR,ZZ,v_t,r_t
      real*8 kappa,C,betaPlus,betaMinus
      real*8 side,xi1Max,xi1Min
      real*8 xi1Plus,xi1Minus,xiStart,xiEnd
      real*8 grav,rho0,complete
      real   completeS,muSingle
      real*8 Phi0,mu,coef,tempbeta,AAA,BBB
      real*8 viPlus,viMinus,xi2Plus,xi2Minus,thetaMax,thetaMin
      real*8 xx,ss,term1,darg1,darg2,tMax,tMin,sumPot
      real   arg1,arg2
      real*8 volMax,volMin,sumVol,Vol0,tempsum
      real*8 xi2One,xi2minOne,volOne,volminOne,volAdd
      real*8 xi(5000),xihalf(4999)
      real*8 pieceIII5,pieceIII3
      integer n,idiag,nzones,iiZone
  201 format(1x,'RR, ZZ, v_t, r_t:',1p4d13.5)
  202 format(1x,'side, kappa, C, betaPlus, betaMinus, xi1Max, xi1Plus,&
     & xi1Minus, xi1Min')
  203 format(1x, 1p9d13.5)
      idiag = 0
      pii = 4.0d0*datan(1.0d0)
      complete = pii/2.0d0
      completeS = complete
      grav = 1.0d0
      rho0 = 1.0d0
      Phi0 = 4.0d0*dsqrt(2.0d0)*grav*rho0*RR**2/3.0d0
      Vol0 = 2.0d0*pii**2*v_t*r_t**2
      side = (v_t**2 - r_t**2)
      kappa = ZZ**2 + RR**2 - side
      C = 1.0d0 + (2.0d0*ZZ/kappa)**2*side
      betaPlus = -(kappa/2.0d0)*(v_t - r_t*dsqrt(C))/side
      betaMinus = -(kappa/2.0d0)*(v_t + r_t*dsqrt(C))/side
</pre>
  </td>
  <td align="left" width="50%">
<ul>
<li>grav --> <math>~G = 1</math></li>
<li>rho0 --> <math>~\rho_0 = 1</math></li>
<li>Phi0 = 4.0d0*dsqrt(2.0d0)*grav*rho0*RR**2/3.0d0 --> <math>~\frac{2^{5 / 2}}{3}\biggl(G\rho_0 R^2\biggr)</math></li>
<li>Vol0 = 2.0d0*pii**2*v_t*r_t**2 &nbsp; &nbsp; <math>~~~~\rightarrow ~~~~ 2\pi^2 \varpi_t r_t^2</math></li>
<li>side = (v_t**2 - r_t**2) &nbsp; &nbsp; <math>~~~~\rightarrow~~~~ (\varpi_t^2 - r_t^2)</math></li>
<li>kappa = ZZ**2 + RR**2 - side &nbsp; &nbsp; <math>~~~~\rightarrow~~~~ \kappa = Z^2 + R^2 - (\varpi_t^2 - r_t^2)</math></li>
<li>C = 1.0d0 + (2.0d0*ZZ/kappa)**2*side &nbsp; &nbsp;  <math>~\rightarrow~~~~ C = 1 + \biggl[\frac{2Z}{Z^2 + R^2 - (\varpi_t^2 - r_t^2)} \biggr]^2(\varpi_t^2 - r_t^2)</math></li>
<li>betaPlus = -(kappa/2.0d0)*(v_t - r_t*dsqrt(C))/side &nbsp; &nbsp;  <math>~\beta_+ = - \frac{\kappa (\varpi_t - r_t \sqrt{C})}{ 2(\varpi_t^2 - r_t^2)}</math></li>
<li>betaMinus = -(kappa/2.0d0)*(v_t + r_t*dsqrt(C))/side &nbsp; &nbsp;  <math>~\beta_- = - \frac{\kappa (\varpi_t + r_t \sqrt{C})}{ 2(\varpi_t^2 - r_t^2)}</math></li>
</ul>
  </td>
</tr>
<tr>
  <td align="left">
<pre>
      if(kappa.le.0.0d0)xi1Max = 1.0d0/dsqrt(1.0d0-(RR/(v_t-betaPlus))**2)
      if(kappa.gt.0.0d0)xi1Min = 1.0d0/dsqrt(1.0d0-(RR/(v_t-betaPlus))**2)
      pieceIII5 = RR**3*xi1Min/(2.0d0*v_t*r_t**2*(xi1Min**2-1.0d0)**1.5)
      if(kappa.le.0.0d0)xi1Min = 1.0d0/dsqrt(1.0d0-(RR/(v_t-betaMinus))**2)
      if(kappa.gt.0.0d0)xi1Max = 1.0d0/dsqrt(1.0d0-(RR/(v_t-betaMinus))**2)
      xi1Plus   = 0.0d0
      xi1Minus  = 0.0d0
      pieceIII3 = 0.0d0
      if(r_t**2.ge.ZZ**2)then
      xi1Plus  = DABS(((v_t + dsqrt(r_t**2 - ZZ**2))**2 + RR**2)/  &
     &           ((v_t + dsqrt(r_t**2 - ZZ**2))**2 - RR**2))
      endif
      if(r_t**2.ge.ZZ**2)then
      xi1Minus = DABS(((v_t - dsqrt(r_t**2 - ZZ**2))**2 + RR**2)/  &
     &           ((v_t - dsqrt(r_t**2 - ZZ**2))**2 - RR**2))
      endif
      if(r_t**2.ge.ZZ**2)then
      pieceIII3 = RR**3*xi1Minus/(2.0d0*v_t*r_t**2*(xi1Minus**2-1.0d0)**1.5)
      endif
!     write(*,201)RR,ZZ,v_t,r_t
!     write(*,202)
!     write(*,203)side,kappa,C,betaPlus,betaMinus,xi1Max,xi1Plus,xi1Minus,xi1Min
</pre>
  </td>
  <td align="left">
<math>~\Rightarrow~~~\mathrm{pieceIII5} = \biggl(\frac{R^3  }{ 2\varpi_t r_t^2}\biggr) \biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_-)^2} \biggr]^{-1 / 2} \biggl\{ \biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_-)^2} \biggr]^{-1 / 2} - 1 \biggr\}^{- 3 / 2}</math>

<p>&nbsp;</p>
<p>&nbsp;</p>
<table border="1" cellpadding="5" align="center">
<tr>
  <td align="center">&nbsp;</td>
  <td align="center"><math>~\xi_1|_\mathrm{min}</math></td>
  <td align="center"><math>~\xi_1|_\mathrm{max}</math></td>
</tr>
<tr>
  <td align="center"><math>~\kappa \le 0</math></td>
  <td align="center"><math>~\biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_-)^2} \biggr]^{-1 / 2}</math></td>
  <td align="center"><math>~\biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_+)^2} \biggr]^{-1 / 2}</math></td>
</tr>
<tr>
  <td align="center"><math>~\kappa > 0</math></td>
  <td align="center"><math>~\biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_+)^2} \biggr]^{-1 / 2}</math></td>
  <td align="center"><math>~\biggl[ 1 - \frac{R^2}{(\varpi_t - \beta_-)^2} \biggr]^{-1 / 2}</math></td>
</tr>
</table>

<p>&nbsp;</p>
<p>&nbsp;</p>
<table border="1" cellpadding="5" align="center">
<tr>
  <td align="center">&nbsp;</td>
  <td align="center"><math>~\xi_1|_-</math></td>
  <td align="center"><math>~\xi_1|_+</math></td>
  <td align="center">pieceIII3</td>
</tr>
<tr>
  <td align="center"><math>~r_t^2 < Z^2</math></td>
  <td align="center"><math>~0</math></td>
  <td align="center"><math>~0</math></td>
  <td align="center"><math>~0</math></td>
</tr>
<tr>
  <td align="center"><math>~r_t^2 \ge Z^2</math></td>
  <td align="center"><math>~~~ \biggl| \frac{ \varpi_t - \sqrt{(r_t^2 - Z^2)^2 + R^2} }{ \varpi_t - \sqrt{ (r_t^2 - Z^2)^2 - R^2 } } \biggr|~~~</math></td>
  <td align="center"><math>~~~\biggl| \frac{ (\varpi_t + \sqrt{ r_t^2 - Z^2})^2 + R^2 }{ ( \varpi_t + \sqrt{r_t^2 - Z^2 } )^2 - R^2} \biggr|~~~</math></td>
  <td align="center"><math>~R^3 \xi_1|_- \biggl[ 2\varpi_t r_t^2 ( \xi^2_1|_- - 1)^{3 / 2}\biggr]^{- 1}</math></td>
</tr>
</table>
  </td>
</tr>
<tr>
  <td align="left">
<pre>
! Begin 1D integration
! Specify sub-region...
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     On 29 May 2016, modified iiZone definition to correspond with
!     online VisTrails chapter discussion:
!     2DStructure/ToroidalCoordinateIntegrationLimits
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! ZONE II
! [iiZone = 1] Green cropped-top sub-volume
      if(iiZone .eq. 1)then
      xi1Start = xi1Minus
      xi1End = xi1Plus
! [iiZone = 3] Lower (orange) Sub-volume that shares edge with 
!              Green cropped-top sub-volume
      else if(iiZone .eq. 3)then
      xi1Start = xi1Minus
      xi1End = xi1Plus
! [iiZone = 4] (Blue) Segment to the left of Lower Sub-volume
      else if(iiZone .eq. 4)then
      xi1Start = xi1Min
      xi1End = xi1Minus
! [iiZone = 2] (Yellow) Segment to the right of Lower Sub-volume
      else if(iiZone .eq. 2)then
      xi1Start = xi1Plus
      xi1End = xi1Max
!
! ZONE III
! [iiZone = 31] (lightblue) Segment of Zone III
      else if(iiZone .eq. 31)then
      xi1End = xi1Max
      xi1Start = xi1Plus
! [iiZone = 32] (darkgreen) Segment of Zone III
      else if(iiZone .eq. 32)then
      xi1End = xi1Plus
      xi1Start = xi1Minus
! [iiZone = 33] (pink=4a) Segment of Zone III
      else if(iiZone .eq. 33)then
      xi1End = xi1Plus
      xi1Start = xi1Min
! [iiZone = 34] (red=4b) Segment of Zone III
      else if(iiZone .eq. 34)then
      xi1End = xi1Minus
      xi1Start = xi1Min
!
! ZONE I
! Full volume
      else
      xi1Start = xi1Min
      xi1End = xi1Max
      end if
! Finished specifying sub-region
</pre>
  </td>
  <td align="left">
Next comments.
  </td>
</tr>
<tr>
  <td align="left">
<pre>
      sumPot = 0.0d0
      sumVol = 0.0d0
      nzones = 1000
      dxi = (xi1End-xi1Start)/dfloat(nzones-1)
      xi(1) = xi1Start
      do n=2,nzones
      xi(n) = xi(n-1)+dxi
      end do
      do n=1,nzones-1
      xihalf(n) = 0.5d0*(xi(n) + xi(n+1))
      end do
! For each value of xi, evaluate these quantities
!!!!Diagnostics!!!!
  303 format(1x,'arg1      = ',1pd15.7,/,&
     &       1x,'arg2      = ',1pd15.7,/,&
     &       1x,'tMax      = ',1pd15.7,/)
  304 format(1x,'arg1      = ',1pd15.7,/,&
     &       1x,'arg2      = ',1pd15.7,/,&
     &       1x,'tMin      = ',1pd15.7,/)
  403 format(1x,'arg1      = ',1pd15.7,/,&
     &       1x,'volMax    = ',1pd15.7,/)
  404 format(1x,'arg1      = ',1pd15.7,/,&
     &       1x,'volMin    = ',1pd15.7,/)
  405 format(1x,'TorusPot  = ',1pd15.7,/,&
     &       1x,'sum       = ',1pd15.7,/,&
     &       1x,'Vol0      = ',1pd15.7,/,&
     &       1x,'sumVol    = ',1pd15.7,/)
! If Z0 less than radius of torus, carry out double integration 
      do n=1,nzones-1
        xx = xihalf(n)
!       xx = xi(n)
        ss = dsqrt(xx**2-1.0d0)
        mu = dsqrt(2.0d0*ss/(ss+xx))
        muSingle = mu
!       coef = dsqrt(xx+1)*ellf(completeS,muSingle)/(ss**4*dsqrt(ss+xx))
        coef = 1.0d0
        tempbeta = v_t/RR - xx/ss
        AAA = (ZZ/RR)**2 + tempbeta**2
        BBB = 2.0d0*v_t*ZZ**2/(RR*kappa) - tempbeta
        term1 = 1.0d0 - AAA*C/BBB**2
        if(term1.lt.0.0d0)term1=0.0d0
        viPlus = kappa*BBB/(2.0d0*RR**2*AAA)*(1.0d0+dsqrt(term1))
        viMinus = kappa*BBB/(2.0d0*RR**2*AAA)*(1.0d0-dsqrt(term1))
        xi2Plus = xx - ss/viPlus
        xi2Minus = xx - ss/viMinus
        thetaMax = dacos(xi2Plus)
        thetaMin = dacos(xi2Minus)
        darg2 = dsqrt(2.0d0/(1.0d0+xx))
        arg2 = darg2
! Determine T(thetaMax)...
!       darg1 = (pii - thetaMax)/2.0d0 
!       arg1 = darg1
!       tMax = dsin(thetaMax)*(5.0d0*xx**2 - 4.0d0*xx*xi2Plus - 1.0d0)&
!    &         /(dsqrt(xx + 1.0d0)*dsqrt((xx-xi2Plus)**3))&
!    &         -4.0d0*xx*elle(arg1,arg2) + (xx-1.0d0)*ellf(arg1,arg2)
!     write(*,310)ss,mu,tempbeta,AAA,BBB,viPlus,viMinus
  310 format(1x,'ss, mu, tempbeta, AAA, BBB, viPlus, viMinus:',&
     &       /,1x,1p7d15.7)
!     write(*,303)arg1,arg2,tMax 
! Determine T(thetaMin)...
!       darg1 = (pii - thetaMin)/2.0d0 
!       arg1 = darg1
!       tMin = dsin(thetaMin)*(5.0d0*xx**2 - 4.0d0*xx*xi2Minus - 1.0d0)&
!    &         /(dsqrt(xx + 1.0d0)*dsqrt((xx-xi2Minus)**3))&
!    &         -4.0d0*xx*elle(arg1,arg2) + (xx-1.0d0)*ellf(arg1,arg2)
!     write(*,304)arg1,arg2,tMin 
!!
!       sumPot = sumPot + dxi*coef*(tMax - tMin)
! TORUS VOLUME DETERMINATION...
        xi2One = 1.0d0
        xi2minOne = -1.0d0
! Determine Vol(thetaMax)...
        volMax = dsqrt(1.0d0-xi2Plus**2)*(4.0d0*xx**2 - 3.0d0*xx*xi2Plus - 1.0d0)&
     &         /(ss**4*(xx-xi2Plus)**2)&
     &         + ((2.0d0*xx**2+1.0d0)/ss**5)*dacos((xx*xi2Plus - 1.0d0)/(xx-xi2Plus)) 
! Determine Vol(thetaMin)...
        volMin = dsqrt(1.0d0-xi2Minus**2)*(4.0d0*xx**2 - 3.0d0*xx*xi2Minus - 1.0d0)&
     &         /(ss**4*(xx-xi2Minus)**2)&
     &         + ((2.0d0*xx**2+1.0d0)/ss**5)*dacos((xx*xi2Minus - 1.0d0)/(xx-xi2Minus)) 
! Determine Vol(One)...
        volOne = 0.0
!       volOne = dsqrt(1.0d0-xi2One**2)*(4.0d0*xx**2 - 3.0d0*xx*xi2One - 1.0d0)&
!    &         /(ss**4*(xx-xi2One)**2)&
!    &         + ((2.0d0*xx**2+1.0d0)/ss**5)*dacos((xx*xi2One - 1.0d0)/(xx-xi2One)) 
! Determine Vol(minusOne)...
        volminOne =  pii*((2.0d0*xx**2+1.0d0)/ss**5) 
!       volminOne = dsqrt(1.0d0-xi2minOne**2)*(4.0d0*xx**2 - 3.0d0*xx*xi2minOne - 1.0d0)&
!    &         /(ss**4*(xx-xi2minOne)**2)&
!    &         + ((2.0d0*xx**2+1.0d0)/ss**5)*dacos((xx*xi2minOne - 1.0d0)/(xx-xi2minOne)) 
! Specify sub-region...
! Full volume or similar segments ...
      if(iiZone.eq.0 .or. iiZone.eq.2 .or. iiZone.eq.4)volAdd = volMax-volMin
! [iiZone = 1] Green cropped-top sub-volume
      if(iiZone .eq. 1)volAdd = volOne - volMax
! [iiZone = 2] Lower (orange) Sub-volume that shares edge with Green cropped-top sub-volume
      if(iiZone .eq. 3)volAdd = volOne-VolMin
! [iiZone = 31] Blue sub-volume of zone III
      if(iiZone.eq.31)volAdd = volMax - volMin
! [iiZone = 32 ] Darkgreen Sub-volumes of zone III
      if(iiZone .eq. 32)volAdd = volOne - volMin 
! [iiZone = 33] Pink Sub-volumes of zone III
      if(iiZone .eq. 33)volAdd = volOne - volMax
! [iiZone = 34] Red Sub-volume zone III
      if(iiZone .eq. 34)volAdd = volMin - volminOne  
        sumVol = sumVol + dxi*volAdd
      enddo
        tempsum = sumVol
        sumVol = pii*RR**3*sumVol/Vol0
!     write(*,405)TorusPot,tempsum,Vol0,sumVol
      return
      end
</pre>
  </td>
  <td align="left">
Next comments.
  </td>
</tr>
</table>
</div>
&nbsp;<br />
<font color="red"><b>Attention!</b></font> &nbsp; As a default, only a portion of the code is displayed in the above table; using your mouse, point anywhere inside the table then scroll down/up to view the entire file.
&nbsp;<br />
&nbsp;<br />

===XML Writer for Visualization===

<div id="XMLwriterForVisualization" style="width: 85%; height: 50em; overflow: auto;">
<table border="1" align="center" width="90%">
<tr>
  <th align="center">Unedited ''XMLwriter01.for''</th>
  <th align="center">Current Comments &amp; Notes</th>
</tr>
<tr>
  <td align="left">
<pre>
      Subroutine XMLwriter01(imax,x,y,cell_scalar)
      real x(50),y(50),z(1)
      real cell_scalar(49,49),point_scalar(50,50)
      integer imax
      integer extentX,extentY,extentZ
      integer ix0,iy0,iz0
      integer norm(50,3)
!     imax=50
      ix0=0
      iy0=0
      iz0=0
      extentX=imax-1
      extentY=imax-1
      extentZ=0
      z(1) = 0.0
! Set normal vector 1D array
      do i=1,imax
        norm(i,1)=0
        norm(i,2)=0
        norm(i,3)=1
      enddo
! Set values of point_scalar array
      do i=1,imax-1
        do j=1,imax-1
        point_scalar(i,j)=cell_scalar(i,j)
        enddo
      enddo
      do i=1,imax-1
        point_scalar(i,imax)=cell_scalar(i,imax-1)
      enddo
      do j=1,imax-1
        point_scalar(imax,j)=cell_scalar(imax-1,j)
      enddo
        point_scalar(imax,imax)=point_scalar(imax-1,imax)
! End point_scalar definition
  201 format('<?xml version="1.0"?>')
  202 format('<VTKFile type="RectilinearGrid" version="0.1" byte_order="LittleEndian">')
  302 format('</VTKFile>')
  203 format(2x,'<RectilinearGrid WholeExtent="',6I3,'">')
  303 format(2x,'</RectilinearGrid>')
  204 format(4x,'<Piece Extent="',6I3,'">')
  304 format(4x,'</Piece>')
  205 format(6x,'<CellData Scalars="cell_scalars" Normals="magnify">')
  305 format(6x,'</CellData>')
  206 format(8x,'<DataArray type="Float32" Name="magnify" NumberOfComponents="3" format="ascii">')
  207 format(8x,'<DataArray type="Float32" Name="cell_scalars" format="ascii">')
  399 format(8x,'</DataArray>')
  208 format(6x,'<PointData Scalars="colorful" Normals="direction">')
  308 format(6x,'</PointData>')
  209 format(8x,'<DataArray type="Float32" Name="colorful" format="ascii">')
  210 format(6x,'<Coordinates>')
  310 format(6x,'</Coordinates>')
  211 format(8x,'<DataArray type="Float32" format="ascii" RangeMin="0" RangeMax="5">')
  212 format(8x,'<DataArray type="Float32" format="ascii">')
  213 format(8x,'<DataArray type="Float32" Name="direction" NumberOfComponents="3" format="ascii">')
  501 format(10f9.5)
  502 format(10f9.5)
  503 format(5x,9(1x,3I2))
  504 format(10f9.5)
  505 format(5x,10(1x,3I2))
!!!!!
!
!  Begin writing out XML tags.
!
!!!!!
      write(*,201)                              !<?xml
      write(*,202)                              !VTKFile
        write(*,203)ix0,extentX,iy0,extentY,iz0,extentZ        !  RectilinearGrid
          write(*,204)ix0,extentX,iy0,extentY,iz0,extentZ      !    Piece
            write(*,205)                        !      CellData
              write(*,207)                      !        DataArray(cell_scalars)
      do j=1,imax-1
      write(*,501)(cell_scalar(i,j),i=1,imax-1)
      enddo
              write(*,399)                      !        /DataArray
              write(*,206)                      !        DataArray(cell_scalars)
      do j=1,imax-1
      write(*,503)(norm(i,1),norm(i,2),norm(i,3),i=1,imax-1)
      enddo
              write(*,399)                      !        /DataArray
            write(*,305)                        !      /CellData
            write(*,208)                        !      PointData
              write(*,209)                      !        DataArray(points)
      write(*,502)((point_scalar(i,j),j=1,imax),i=1,imax)
              write(*,399)                      !        /DataArray
              write(*,213)                      !        DataArray(cell_scalars)
      do j=1,imax
      write(*,505)(norm(i,1),norm(i,2),norm(i,3),i=1,imax)
      enddo
              write(*,399)                      !        /DataArray
            write(*,308)                        !      /PointData
            write(*,210)                        !      Coordinates
              write(*,212)                      !        DataArray(x-direction)
        write(*,504)(x(i),i=1,imax)
              write(*,399)                      !        /DataArray
              write(*,212)                      !        DataArray(y-direction)
        write(*,504)(y(i),i=1,imax)
              write(*,399)                      !        /DataArray
              write(*,212)                      !        DataArray(z-direction)
        write(*,504)z(1)
              write(*,399)                      !        /DataArray
            write(*,310)                        !      /Coordinates
          write(*,304)                          !    /Piece
        write(*,303)                            !  /RectilinearGrid
      write(*,302)                              !/VTKFile
      return
      end
</pre>
  </td>
  <td align="left" width="50%">
<ul>
<li>dvarpi_t: &nbsp; <math>~\varpi_t = \tfrac{3}{4}</math></li>
</ul>
  </td>
</tr>
</table>
</div>
&nbsp;<br />
<font color="red"><b>Attention!</b></font> &nbsp; As a default, only a portion of the code is displayed in the above table; using your mouse, point anywhere inside the table then scroll down/up to view the entire file.
&nbsp;<br />
&nbsp;<br />