Editing 2DStructure/ToroidalCoordinateIntegrationLimits (section)

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