Editing 2DStructure/ToroidalCoordinateIntegrationLimits (section)

===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 />