c plot Modular Coil Section Measurement Data
c plotmeas.f  - plot After Initial Winding Measurements as Coil Cross Section
c               and Target Cross Section based on Mike Z's Area Preserving Algorithm
c plotmeas2.f - modified to calculate target based on Minimum Change Algorithm
c plotmeas3.f - modified to calculate target based on Minimum Change Algorithm
c               and also target inputted Current Center
c plotmeas4.f - revert to Area Preserving Algorithm
c               and also target inputted Current Center
c                Add option iadjust
c           =1: do Area Preserve, =2: Do minimum change, =3: Initial Settings/Targets=0: Plot All
c plotmeas5.f - Revise code to optional calculate Initial Settings/Targets
c               based on inputted coil curvatures and torsion
c          Clean up redundant coding
c          Add flags to control execution
c*******************************************************************************
      module measurements
c*******************************************************************************
      implicit real*8(a-h,o-z)
      real*8, dimension(9) :: x0, y0 ! Winding Form
      real*8, dimension(100) :: xA1, yA1, xB1, yB1 ! Windings (undeformed)
      real*8, dimension(100) :: xA2, yA2, xB2, yB2 ! Windings (deformed)
      real*8, dimension(100) :: xA3, yA3, xB3, yB3 ! Windings (adjusted target Area Preserving)
      real*8, dimension(100) :: xA4, yA4, xB4, yB4 ! Windings (adjusted target Min Change)
      real*8, dimension(100) :: xA5, yA5, xB5, yB5 ! Windings (Initial Settings/Targets)
      real*8, dimension(:), allocatable  :: dxc, dyc ! Current Center Target
      real*8, dimension(:,:), allocatable :: baseA, teeA, baseB, teeB  ! Initial Measurements of MCWF
      real*8, dimension(:,:), allocatable :: topA, sideA, topB, sideB  ! Initial Meas after Winding
      real*8, dimension(:,:), allocatable ::topA3, sideA3, topB3, sideB3 ! Windings (adjusted target Area Preserving)
      real*8, dimension(:,:), allocatable ::topA4, sideA4, topB4, sideB4 ! Windings (adjusted target Min Change)
      real*8, dimension(:,:), allocatable ::topA5, sideA5, topB5, sideB5 ! Windings (Initial Settings/Targets)
      real*8, dimension(:), allocatable :: Curv_tot, Curv_InPl, Curv_OOP
     +                                    , Tors
      character*80 title, labx, laby,  subtitle, labz
      character*32 baseAA, teeAA, topAA, sideAA
      character*32 baseBB, teeBB, topBB, sideBB
      character*32  curcen, initial, coilcurv
      character*4 name(100), cdum
      character*16 lab
      real*8, dimension(2) :: x, y
      logical :: LreadCurv, LreadCurCen, LreadMCWF, LreadWind
      integer :: n1, n2, n3, n4, nsect, nw, nh, nclamps
      namelist/pltdat/ title, subtitle, labx, laby, labz,
     + xmin, xmax, ymin, ymax,  h, w, t, b, scale, dr,
     +nspc, nw, nh, nclamps, iadjust,
     + baseAA, teeAA, topAA, sideAA, baseBB, teeBB, topBB, sideBB,
     +  curcen, initial, coilcurv,
     +  LreadCurv, LreadCurCen, LreadMCWF, LreadWind
      end module measurements
c
c*******************************************************************************
      program main
c*******************************************************************************
      use measurements
      implicit real*8(a-h,o-z)
c
c Read in Namelist data
      read(*,pltdat)
      nsect=nclamps*nspc ! nspc == #sect per clamp -> USE 1
c
c Allocate Arrays
c
      allocate(dxc(nsect), dyc(nsect))
      allocate(Curv_tot(nsect), Curv_InPl(nsect))
      allocate(Curv_OOP(nsect), Tors(nsect))
c
      allocate(baseA(nw,nsect))
      allocate(baseB(nw,nsect))
      allocate(topA(nw,nsect))
      allocate(topB(nw,nsect))
c
      allocate(teeA(nh,nsect))
      allocate(teeB(nh,nsect))
      allocate(sideA(nh,nsect))
      allocate(sideB(nh,nsect))
c
      allocate(topA3(nw,nsect))
      allocate(topB3(nw,nsect))
      allocate(sideA3(nh,nsect))
      allocate(sideB3(nh,nsect))
c
      allocate(topA4(nw,nsect))
      allocate(topB4(nw,nsect))
      allocate(sideA4(nh,nsect))
      allocate(sideB4(nh,nsect))
c
      allocate(topA5(nw,nsect))
      allocate(topB5(nw,nsect))
      allocate(sideA5(nh,nsect))
      allocate(sideB5(nh,nsect))
c
c  Read in Current Center Target if requested
c
      if(LreadCurCen)then
        open(unit=11,file=curcen,status="old")
        read(11,*)
        read(11,*)
        read(11,*)
        do i=1,nsect
        read(11,*) idum, dxc(i), dyc(i)
        enddo
        close(11)
      else
        dxc = 0.
        dyc = 0.
      endif
c
c  Read in Coil Curvatures and Torsion and Calcualte Initial Settings
c  or Read in precalculated Initial Settings
c
      if(LreadCurv)then  ! Read in Coil Curvatures and Torsion
        open(unit=11,file=coilcurv,status="old")
        read(11,*)
        do i=1,nsect
        read(11,*) cdum, idum2,
     +          Curv_tot(i), Curv_InPl(i), Curv_OOP(i), Tors(i)
        enddo
        close(11)
      else  ! Read in Initial Side Clamp Settings and Targets
        open(unit=5,file=initial,status="old")
        read(5,*)
        do isect=1,nsect
        read(5,*) w1, h1, w2, h2
        sideA5(:,isect)=-w1
        sideB5(:,isect)= w2
        topA5(:,isect) = h1
        topB5(:,isect) = h2
        enddo
        close(5)
      endif
c
c
c FIX to Fill in clamp Names
c
      do i=1,95
      ii=i
      if(i.gt.14)ii=i+1
      write(name(i),"(i2)") ii
      enddo
      name(14)=trim(name(14))//"a"
c
c Read in romer arm data for Winding Form Measurements
c
      if(LreadMCWF)then
        open(unit=5,file=baseAA,status="old")
        do isect=1,nsect
        read(5,*) idum,( baseA(i,isect), i=1,nw)
        enddo
        close(5)
c
        open(unit=5,file=baseBB,status="old")
        do isect=1,nsect
        read(5,*) idum,( baseB(i,isect), i=1,nw)
        enddo
        close(5)
c
        open(unit=5,file=teeAA,status="old")
        do isect=1,nsect
        read(5,*) idum,( teeA(i,isect),i=1,nh)
        enddo
        close(5)
c
        open(unit=5,file=teeBB,status="old")
        do isect=1,nsect
        read(5,*) idum,( teeB(i,isect),i=1,nh)
        enddo
        close(5)
      endif
c
c
c Read in romer arm data for Initial Measurements after Winding
c
      if(LreadWind)then
        open(unit=5,file=topAA,status="old")
        do isect=1,nsect
        read(5,*) idum,( topA(i,isect), i=1,nw)
        enddo
        close(5)
c
        open(unit=5,file=topBB,status="old")
        do isect=1,nsect
        read(5,*) idum,( topB(i,isect), i=1,nw)
        enddo
        close(5)
c
        open(unit=5,file=sideAA,status="old")
        do isect=1,nsect
        read(5,*) idum,( sideA(i,isect),i=1,nh)
        enddo
        close(5)
c
        open(unit=5,file=sideBB,status="old")
        do isect=1,nsect
        read(5,*) idum,( sideB(i,isect),i=1,nh)
        enddo
        close(5)
      endif
c
c Calculate Initial Side Clamp Settings and Top Clamp Targets if data available
c
      if(LreadMCWF .and. LreadCurv) then
          call InitialSet
      endif
c
c Calculate New Targets based on Initial Measurements
c
!      if(iadjust.eq.1)then
          call readjust3
!      elseif(iadjust.eq.2)then
          call readjust4
!      endif
c
c Scale Deformations for Measured Data and Targets for Plots
c
      baseA=baseA*scale
      baseB=baseB*scale
      topA=topA*scale
      topB=topB*scale
      teeA=teeA*scale
      teeB=teeB*scale
      sideA=sideA*scale
      sideB=sideB*scale
c
      topA3=topA3*scale
      topB3=topB3*scale
      sideA3=sideA3*scale
      sideB3=sideB3*scale
c
      topA4=topA4*scale
      topB4=topB4*scale
      sideA4=sideA4*scale
      sideB4=sideB4*scale
c
      topA5=topA5*scale
      topB5=topB5*scale
      sideA5=sideA5*scale
      sideB5=sideB5*scale
c
c Initialize Winding Form and Undeformed Winding Shape
c
c
c Order points as base, tee, top, side
c
      n1= 1+nw+1
      n2=n1+nh+1
      n3=n2+nw+1
      n4=n3+nh+1
c
      xx1=-w-t/2
      xx2=-t/2
      xx3=t/2
      xx4=t/2+w
      yy1=-h/2-b
      yy2=-h/2
      yy3=h/2
      x0=(/xx1,xx1,xx2,xx2,xx3,xx3,xx4,xx4,xx1/)
      y0=(/yy1,yy2,yy2,yy3,yy3,yy2,yy2,yy1,yy1/)
c Order points as base, tee, top, side
      dx=w/nw
      dy=h/nh
c
      xA1(1)=xx1
      xA1(2:n1-1)=(/ (xx1+dx/2+i*dx, i=0,nw-1) /)
      xA1(n1:n2)=xx2
      xA1(n2+1:n3-1)=(/ (xx2-dx/2-i*dx, i=0,nw-1) /)
      xA1(n3:n4)=xx1
c
      yA1(1:n1)=yy2
      yA1(n1+1:n2-1)=(/ (yy2+dy/2+i*dy, i=0,nh-1) /)
      yA1(n2:n3)= yy3
      yA1(n3+1:n4-1)=(/ (yy3-dy/2-i*dy, i=0,nh-1) /)
      yA1(n4)=yy2
c
      xB1(:)=-xA1(:)
      yB1(:)=+yA1(:)
c
      call plt0
      isect=0
      do i=1,nclamps
      do j=1,1 ! nspc
!      isect=isect+1
      isect=1+(i-1)*nspc
c
c----------------------------------------------------------------------------
c
c Add Measured Cross section
c
      call Fill_Xsect(topA,sideA,topB,sideB,xA2,yA2,xB2,yB2,isect)
c----------------------------------------------------------------------------
c
c Add Adjusting Target Section - Area Preserving
c
      call Fill_Xsect(topA3,sideA3,topB3,sideB3,xA3,yA3,xB3,yB3,isect)
c----------------------------------------------------------------------------
c
c Add Adjusting Target Section - Minimum Change
c
      call Fill_Xsect(topA4,sideA4,topB4,sideB4,xA4,yA4,xB4,yB4,isect)
c----------------------------------------------------------------------------
c
c Add Adjusting Target Section - Initial targets
c
      call Fill_Xsect(topA5,sideA5,topB5,sideB5,xA5,yA5,xB5,yB5,isect)
c----------------------------------------------------------------------------
c
      write(subtitle,"('Clamp ',a4,' Section ',i1,'$')") name(i),j
      call plt1(xmin,xmax,ymin,ymax,title,subtitle,labx,laby,"$")
      call plt2(x0, y0, 9, 1, "$")     ! Winding Form
      call plt2(xA1, yA1, n4, 1, "$")  ! Windings (undeformed)
      call plt2(xB1, yB1, n4, 1, "$")
      call plchhq2(-1.d0, 0.d0,"A$",.02d0,0.d0,0)
      call plchhq2(1.d0, 0.d0,"B$",.02d0,0.d0,0)
      call plt2(xA2, yA2, n4, 2, "$")  ! Windings (deformed)
      call plt2(xB2, yB2, n4, 2, "$")
      if(iadjust.eq.1)then
        call plt2(xA3, yA3, n4, 4, "$")  ! Windings (adjusted target Area Preserving)
        call plt2(xB3, yB3, n4, 4, "$")
        write(lab,"(f6.3,'$')") sideA3(1,isect)/scale
        call plchhq2(-3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") sideB3(1,isect)/scale
        call plchhq2(+3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topA3(1,isect)/scale
        call plchhq2(-1.2d0,1.5d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topB3(1,isect)/scale
        call plchhq2(+1.2d0,1.5d0,lab,.02d0,0.d0,0)
      elseif(iadjust.eq.2)then
        call plt2(xA4, yA4, n4, 3, "$")  ! Windings (adjusted target Min Change)
        call plt2(xB4, yB4, n4, 3, "$")
        write(lab,"(f6.3,'$')") sideA4(1,isect)/scale
        call plchhq2(-3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") sideB4(1,isect)/scale
        call plchhq2(+3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topA4(1,isect)/scale
        call plchhq2(-1.2d0,1.5d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topB4(1,isect)/scale
        call plchhq2(+1.2d0,1.5d0,lab,.02d0,0.d0,0)
      elseif(iadjust.eq.3)then
        call plt2(xA5, yA5, n4, 6, "$")  ! Windings (Initial settings/targets)
        call plt2(xB5, yB5, n4, 6, "$")
        write(lab,"(f6.3,'$')") sideA5(1,isect)/scale
        call plchhq2(-3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") sideB5(1,isect)/scale
        call plchhq2(+3.d0, 0.d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topA5(1,isect)/scale
        call plchhq2(-1.2d0,1.5d0,lab,.02d0,0.d0,0)
        write(lab,"(f6.3,'$')") topB5(1,isect)/scale
        call plchhq2(+1.2d0,1.5d0,lab,.02d0,0.d0,0)
      else ! do all
        call plt2(xA3, yA3, n4, 4, "$")  ! Windings (adjusted target Area Preserving)
        call plt2(xB3, yB3, n4, 4, "$")
        call plt2(xA4, yA4, n4, 3, "$")  ! Windings (adjusted target Min Change)
        call plt2(xB4, yB4, n4, 3, "$")
        call plt2(xA5, yA5, n4, 6, "$")  ! Windings (Initial settings/targets)
        call plt2(xB5, yB5, n4, 6, "$")
      endif
c
c add current center
c
      dxs=dxc(i)*scale
      dys=dyc(i)*scale
      x=(/0., 0./)
      y=(/-.2, +.2/)
      call plt2(x, y, 2, 1, "$")
      x=(/-.2, .2/)
      y=(/0., 0./)
      call plt2(x, y, 2, 1, "$")
      x=(/dxs,dxs/)
      y=(/dys-.2,dys+.2/)
      call plt2(x, y, 2, 4, "$")
      x=(/dxs-.2,dxs+.2/)
      y=(/dys,dys/)
      call plt2(x, y, 2, 4, "$")
      call pltadv
      enddo
      enddo
c
      call pltend
      stop
      end program main
c*******************************************************************************
      subroutine readjust3
      use measurements
      implicit real*8(a-h,o-z)
c  reajust side and top clamps using measured data
c  implement Mike Z's area preserving algorithm
c
      h1min = 1.e9
      h2min = 1.e9
      w1min = 1.e9
      w2min = 1.e9
      h1max = -1.e9
      h2max = -1.e9
      w1max = -1.e9
      w2max = -1.e9
c
      open(unit=8,file="meas.sum",status="replace")
      write(*,"(' Area Preserving Algorithm')")
      write(*,"(' Clamp#  Side A  Side B  Top A   Top B')")
      write(8,"(' Clamp#  x1      y1      x2      y2      w1      h1
     +  w2      h2      s1      s2')")
      do i=1,nclamps
      isect=(i-1)*nspc + 1
c
c use average of measurements along each face
c
      x1 = -sum( teeA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
      x2 =  sum( teeB(1:nh,isect))/nh
      w1 = -sum(sideA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
      w2 =  sum(sideB(1:nh,isect))/nh
      y1 =  sum(baseA(1:nw,isect))/nw
      y2 =  sum(baseB(1:nw,isect))/nw
      h1 =  sum( topA(1:nw,isect))/nw
      h2 =  sum( topB(1:nw,isect))/nw
c
c Mike Z's Algorithm (Solution to 4 equations, 4 unknowns)
c   dyc = (y1 + h1 + y2 + h2)/4  ! preserve current center
c   dxc = (x1 + w1 + x2 + w2)/4  ! preserve current center
c   S1 = H*(x1-w1) + W*(h1-y1)  ! preserve area
c   S2 = H*(w2-x2) + W*(h2-y2)  ! preserve area
c
c
      s1 = h*(x1-w1) + w*(h1-y1)
      s2 = h*(w2-x2) + w*(h2-y2)
c
      h1_new = -h/w*(x1+x2) - y2 - (s2-s1)/2./w
     +         + 2.*dyc(i) + 2.*dxc(i)*h/w
      h2_new =  h/w*(x1+x2) - y1 + (s2-s1)/2./w
     +         + 2.*dyc(i) - 2.*dxc(i)*h/w
      w1_new = -w/h*(y1+y2) - x2 - (s2+s1)/2./h
     +         + 2.*dxc(i) + 2.*dyc(i)*w/h
      w2_new =  w/h*(y1+y2) - x1 + (s2+s1)/2./h
     +         + 2.*dxc(i) - 2.*dyc(i)*w/h
c
      sideA3(:,isect)=-w1_new
      sideB3(:,isect)= w2_new
      topA3(:,isect) = h1_new
      topB3(:,isect) = h2_new
c
      write(8,"(a4,1x,10f8.3)") name(i),x1,y1,x2,y2,w1,h1,w2,h2,s1,s2
      write(*,"(a4,1x,4f8.3)") name(i), -w1_new, w2_new, h1_new, h2_new
c
      h1min=min(h1min,h1_new)
      h2min=min(h2min,h2_new)
      w1min=min(w1min,w1_new)
      w2min=min(w2min,w2_new)
      h1max=max(h1max,h1_new)
      h2max=max(h2max,h2_new)
      w1max=max(w1max,w1_new)
      w2max=max(w2max,w2_new)
      enddo
c
      write(*,"('Min  ',4f8.3)")  -w1max, w2min, h1min, h2min
      write(*,"('Max  ',4f8.3)")  -w1min, w2max, h1max, h2max
      return
      end
c*******************************************************************************
C      subroutine readjust4
C      use measurements
C      implicit real*8(a-h,o-z)
Cc  reajust side and top clamps using measured data
Cc  USE MINIMUM ADJUSTMENT FROM EXISTING STATE
Cc   ie calculate equal adjustments (dw=dw1=dw2, dh=dh1=dh2) such that
Cc     (h1+dh)+(h2+dh)+y1+y2 = 0.
Cc     (w1+dw)+(w2+dw)+x1+x2 = 0.
Cc
C      h1min = 1.e9
C      h2min = 1.e9
C      w1min = 1.e9
C      w2min = 1.e9
C      h1max = -1.e9
C      h2max = -1.e9
C      w1max = -1.e9
C      w2max = -1.e9
Cc
C      open(unit=8,file="meas.sum",status="replace")
C      write(*,"(' Clamp#  Side A  Side B  Top A   Top B')")
C      write(8,"(' Clamp#  x1      y1      x2      y2      w1      h1
C     +  w2      h2      s1      s2')")
C      do i=1,nclamps
C      isect=(i-1)*nspc + 1
Cc
Cc use average of measurements along each face
Cc
C      x1 = -sum( teeA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
C      x2 =  sum( teeB(1:nh,isect))/nh
C      w1 = -sum(sideA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
C      w2 =  sum(sideB(1:nh,isect))/nh
C      y1 =  sum(baseA(1:nw,isect))/nw
C      y2 =  sum(baseB(1:nw,isect))/nw
C      h1 =  sum( topA(1:nw,isect))/nw
C      h2 =  sum( topB(1:nw,isect))/nw
Cc
Cc DO MINIMUM ADJUSTMENT FROM EXISTING STATE
Cc
Cc
C      s1 = h*(x1-w1) + w*(h1-y1)
C      s2 = h*(w2-x2) + w*(h2-y2)
Cc
Cc      h1_new = -h/w*(x1+x2) - y2 - (s2-s1)/2./w
Cc      h2_new =  h/w*(x1+x2) - y1 + (s2-s1)/2./w
Cc      w1_new = -w/h*(y1+y2) - x2 - (s2+s1)/2./h
Cc      w2_new =  w/h*(y1+y2) - x1 + (s2+s1)/2./h
Cc
C      dw=-(w1+w2+x1+x2)/2.
C      dh=-(h1+h2+y1+y2)/2.
Cc
C      h1_new = h1+dh
C      h2_new = h2+dh
C      w1_new = w1+dw
C      w2_new = w2+dw
Cc
C      sideA4(:,isect)=-w1_new
C      sideB4(:,isect)= w2_new
C      topA4(:,isect) = h1_new
C      topB4(:,isect) = h2_new
Cc
C      write(8,"(a4,1x,10f8.3)") name(i),x1,y1,x2,y2,w1,h1,w2,h2,s1,s2
C      write(*,"(a4,1x,4f8.3)") name(i), -w1_new, w2_new, h1_new, h2_new
Cc
C      h1min=min(h1min,h1_new)
C      h2min=min(h2min,h2_new)
C      w1min=min(w1min,w1_new)
C      w2min=min(w2min,w2_new)
C      h1max=max(h1max,h1_new)
C      h2max=max(h2max,h2_new)
C      w1max=max(w1max,w1_new)
C      w2max=max(w2max,w2_new)
C      enddo
Cc
C      write(*,"('Min  ',4f8.3)")  -w1max, w2min, h1min, h2min
C      write(*,"('Max  ',4f8.3)")  -w1min, w2max, h1max, h2max
C      return
C      end
c*******************************************************************************
      subroutine readjust4
      use measurements
      implicit real*8(a-h,o-z)
c  reajust side and top clamps using measured data
c  USE MINIMUM ADJUSTMENT FROM EXISTING STATE to TARGET given Current Center Shift
c   ie calculate equal adjustments (dw=dw1=dw2, dh=dh1=dh2) such that
c     (h1+dh)+(h2+dh)+y1+y2 = 4*dyc
c     (w1+dw)+(w2+dw)+x1+x2 = 4*dxc
c
      h1min = 1.e9
      h2min = 1.e9
      w1min = 1.e9
      w2min = 1.e9
      h1max = -1.e9
      h2max = -1.e9
      w1max = -1.e9
      w2max = -1.e9
c
      open(unit=8,file="meas.sum",status="replace")
      write(*,"(' Minimum Change Algorithm')")
      write(*,"(' Clamp#  Side A  Side B  Top A   Top B')")
      write(8,"(' Clamp#  x1      y1      x2      y2      w1      h1
     +  w2      h2      s1      s2')")
      do i=1,nclamps
      isect=(i-1)*nspc + 1
c
c use average of measurements along each face
c
      x1 = -sum( teeA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
      x2 =  sum( teeB(1:nh,isect))/nh
      w1 = -sum(sideA(1:nh,isect))/nh ! RomerArm measures relative to outward surface normal
      w2 =  sum(sideB(1:nh,isect))/nh
      y1 =  sum(baseA(1:nw,isect))/nw
      y2 =  sum(baseB(1:nw,isect))/nw
      h1 =  sum( topA(1:nw,isect))/nw
      h2 =  sum( topB(1:nw,isect))/nw
c
c DO MINIMUM ADJUSTMENT FROM EXISTING STATE
c
c
      s1 = h*(x1-w1) + w*(h1-y1)
      s2 = h*(w2-x2) + w*(h2-y2)
c
c      h1_new = -h/w*(x1+x2) - y2 - (s2-s1)/2./w
c      h2_new =  h/w*(x1+x2) - y1 + (s2-s1)/2./w
c      w1_new = -w/h*(y1+y2) - x2 - (s2+s1)/2./h
c      w2_new =  w/h*(y1+y2) - x1 + (s2+s1)/2./h
c
      dw=-(w1+w2+x1+x2)/2. + 2.*dxc(i)
      dh=-(h1+h2+y1+y2)/2. + 2.*dyc(i)
c
      h1_new = h1+dh
      h2_new = h2+dh
      w1_new = w1+dw
      w2_new = w2+dw
c
      sideA4(:,isect)=-w1_new
      sideB4(:,isect)= w2_new
      topA4(:,isect) = h1_new
      topB4(:,isect) = h2_new
c
      write(8,"(a4,1x,12f8.3)") name(i),x1,y1,x2,y2,w1,h1,w2,h2,s1,s2,
     +   dxc(i),dyc(i)
      write(*,"(a4,1x,4f8.3)") name(i), -w1_new, w2_new, h1_new, h2_new
c
      h1min=min(h1min,h1_new)
      h2min=min(h2min,h2_new)
      w1min=min(w1min,w1_new)
      w2min=min(w2min,w2_new)
      h1max=max(h1max,h1_new)
      h2max=max(h2max,h2_new)
      w1max=max(w1max,w1_new)
      w2max=max(w2max,w2_new)
      enddo
c
      write(*,"('Min  ',4f8.3)")  -w1max, w2min, h1min, h2min
      write(*,"('Max  ',4f8.3)")  -w1min, w2max, h1max, h2max
      return
      end
c*******************************************************************************
c*******************************************************************************
      module pltmod
c*******************************************************************************
      implicit real*8(a-h,o-z)
      real*4 data_x(2),data_y(2), x0,x1, xp1,xp2, yp1,yp2
      real*4 xmin4,xmax4,ymin4,ymax4
      real*4 vxmi4,vxma4,vymi4,vyma4
      real*4, allocatable, dimension(:) :: x4, y4
      integer iwkid
      logical legend
      save
       character*72 text   ! internal22
      data nframe/0/,ncurve/0/,ntot/0/
      end module pltmod
c*******************************************************************************
      subroutine plt0(name)
c*******************************************************************************
c open graphics system
      call gopks(2,idum)
c
c open device
      iwkid=1
      iconid=2 ! fortran unit number
      iwtype =1 ! workstation type - 1=create cgm binary file gmeta
      call gopwk(iwkid,iconid,iwtype)
c
c activate device
      call gacwk(iwkid)
      call dfclrs(iwkid)
c
      open(unit=25,file="pltdat",status="replace")
      return
      end
c
c*******************************************************************************
      subroutine plt1(xmin,xmax,ymin,ymax,title,subtitle,labx,laby,labz)
c*******************************************************************************
c step 3: define picture (plots)
      use pltmod
      implicit real*8(a-h,o-z)
      character*(*) title,subtitle,labx,laby,labz
c
c set color index to black (1)
      call gsplci(1)
      call gstxci(1)
      xmin4=xmin
      xmax4=xmax
      ymin4=ymin
      ymax4=ymax
      vxmi4=.10d0
      vymi4=.10d0
      vxma4=.8d0
      vyma4=.8d0
      if(labz(1:1).eq.'$')vxma4=.90d0
      if(labz(1:1).eq.'$')vyma4=.90d0
      legend=.true.
      if(labz(1:1).eq.'$')legend=.false.
c set coordinate system
      itype=1 ! define transformation type 1=linlin
      call set(vxmi4,vxma4,vymi4,vyma4,xmin4,xmax4,ymin4,ymax4,itype)
c
c
c set world coordinates equate to unity for text
      call getset(vxmi4,vxma4,vymi4,vyma4,xmin4,xmax4,ymin4,ymax4,itype)
c
c set ticks only
      x1=1.d0
      call agsetr('BACKGROUND.',x1) ! sete perimeter background
C      call agsetr('AXIS/TICKS/MINOR/SPACING.',0.d0)
c plot grid ( the hard way )
      data_x(1)=xmin4
      data_x(2)=xmax4
      data_y(1)=ymin4
      data_y(2)=ymax4
      call agseti('frame.',2)
      call agseti('SET.',-4)
      call agsetc('LABEL/NAME.','B')
      call agseti('LABEL/SUPPRESSION FLAG.',1)
      call agsetc('LABEL/NAME.','L')
      call agseti('LABEL/SUPPRESSION FLAG.',1)
      call ezmxy(data_x,data_y,1,1,1,' ')
c
c plot titles and labels
      x0=0.d0
      x1=1.d0
      call set(x0,x1,x0,x1,x0,x1,x0,x1,1)
c     call plchhq2(xpos,ypos,chrs,size,ang,cntr)
      call plchhq2(.5d0,.95d0,title,.02d0,0.d0,0)
      if(subtitle(1:1).ne.'$')then
         call plchhq2(.5d0,.87d0,subtitle,.02d0,0.d0,0)
      endif
c      call plchhq2(.5d0,0.92d0,stamp,.01d0,0.d0,0)
      call plchhq2(.5d0,.02d0,labx,.0154d0,0.d0,0)
      call plchhq2(.01d0,.5d0,laby,.0154d0,90.d0,0)
c
c set grid parameters
c      call agsetc("name","char_value")
c      call agseti("name",int_value)
c      call agsetr("name",real_value)
C      call gsgrxm(0)
C      call gsgrym(0)
c
      if(labz(1:1).ne.'$')then
      vymax=vyma4
      call plchhq2(.82d0,vymax,labz,.008d0,0.d0,0)
      endif
c
      write(25,*) ! blank line will signify start of new frame
      write(25,"(a80)") title
      write(25,"(a80)") subtitle
      write(25,"(a80)") labx
      write(25,"(a80)") laby
      write(25,"(a80)") labz
      call set(vxmi4,vxma4,vymi4,vyma4,xmin4,xmax4,ymin4,ymax4,itype)
      return
      end
c*******************************************************************************
      subroutine plt2(x,y,n,iz,lab)
c*******************************************************************************
      use pltmod
      implicit real*8(a-h,o-z)
      dimension x(n), y(n)
      character*(*) lab
c
      allocate(x4(n),y4(n))
c
      do i=1,n
      x4(i)=x(i)
      y4(i)=y(i)
      enddo
c
c set line style
      call gsln(1) ! good only for iz=1.d0.4
c set line size
      call gslwsc(x1) ! default
c  set color index
      call gsplci(iz)
      call gstxci(iz)
c plot line with n points
      call dpcurv(x4,y4,n)
      write(25,"(i5,a8)") n,lab
      write(25,"(1p2e15.6)") (x4(i),y4(i),i=1,n)
c
      if(legend)then
      call getset(vxmi4,vxma4,vymi4,vyma4,xmin4,xmax4,ymin4,ymax4,itype)
          yp1=vyma4 - .03d0* iz
          call set(x0,x1,x0,x1,x0,x1,x0,x1,1)
          data_x(1)=.82d0
          data_x(2)=.90d0
          data_y(1)=yp1
          data_y(2)=yp1
c          call dpline(xp1,xp2,yp1,yp1)
           call dpcurv(data_x,data_y,2)
          yp=yp1+.01d0
          call plchhq2(.9d0,yp,lab,.008d0,0.d0,0)
      call set(vxmi4,vxma4,vymi4,vyma4,xmin4,xmax4,ymin4,ymax4,itype)
      endif
      deallocate(x4,y4)
      return
      end
c
c*******************************************************************************
      subroutine pltadv
c*******************************************************************************
c
c step 4: display pictures (advance frame for new plots)
      use pltmod
      implicit real*8(a-h,o-z)
      nframe=nframe+1
      ntot=ntot+ncurve
      ncurve=0
      write(text,199) nframe
 199  format('figure ',i2,'$')
CAWB      call getset(vxmi4,vxma4,vymi4,vyma4,xmin,xmax,ymin,ymax,itype)
CAWB          call set(0.d0,1.d0,0.d0,1.d0,0.d0,1.d0,0.d0,1.d0,1)
CAWB      call plchhq2(.82d0,.02d0,text,.0154d0,0.d0,-1)
CAWB      call set(vxmi4,vxma4,vymi4,vyma4,xmin,xmax,ymin,ymax,itype)
c advance frame
      call frame
c
      return
      end
c
c*******************************************************************************
      subroutine pltend
c*******************************************************************************
c step 5: terminate library
      use pltmod
      implicit real*8(a-h,o-z)
c
c deactivate device
      iwkid=1
      call gdawk(iwkid)
c
c close device
      call gclwk(iwkid)
c
c close graphics system
      call gclks
c
      return
      end
c*******************************************************************************
      subroutine plchhq2(xpos,ypos,chrs,size,ang,icntr)
c*******************************************************************************
      implicit real*8(a-h,o-z)
      real*4 xpos4,ypos4,size4,ang4
      character*(*) chrs
c find length of string
      do i=1,80
      if(chrs(i:i).eq."$")goto 10
      enddo
      i=81
 10   continue
      if(i.gt.1)then
        xpos4=xpos
        ypos4=ypos
        size4=size
        ang4=ang
        call plchhq(xpos4,ypos4,chrs(1:i-1),size4,ang4,icntr)
      endif
      return
      end
c*******************************************************************************
      SUBROUTINE DFCLRS(IWKID)
c*******************************************************************************
C
C Define a set of RGB color triples for colors 1 through 15 on
C workstation 1.d0
C
        real*4 RGBV(3,15), one
C
C Define the RGB color triples needed below.
C
        DATA RGBV / 0.00 , 0.00 , 0.00 ,  ! Black
     +              1.00 , 0.00 , 0.00 ,  ! Red
     +              0.00 , 1.00 , 0.00 ,  ! Green
     +              0.00 , 0.00 , 1.00 ,  ! Blue
     +              1.00 , 1.00 , 0.00 ,  ! Yellow
     +              1.00 , 0.00 , 1.00 ,  ! Magenta
     +              0.00 , 1.00 , 1.00 ,  ! Cyan
     +              1.00 , 0.50 , 0.00 ,  ! Orange (Red-Yellow)
     +              0.00 , 1.00 , 0.50 ,  ! Green-Cyan
     +              1.00 , 0.00 , 0.50 ,  ! Red-Magenta
     +              1.00 , 0.75 , 0.00 ,  !
     +              1.00 , 0.38 , 0.38 ,  !
     +              1.00 , 0.00 , 0.00 ,  !
     +              0.70 , 1.00 , 0.00 ,  !
     +              1.00 , 0.00 , 0.38 /  !
      data one/1./
C
C Define 16 different color indices, for indices 0 through 15.  The
C color corresponding to index 0 is black and the color corresponding
C to index 1 is white.
C
        CALL GSCR (IWKID,0,one,one,one)
C
        DO 101 I=1,15
          CALL GSCR (IWKID,I,RGBV(1,I),RGBV(2,I),RGBV(3,I))
  101   CONTINUE
C
C Done.
C
        RETURN
C
      END
c*******************************************************************************
      subroutine Fill_Xsect(topA,sideA,topB,sideB,xA2,yA2,xB2,yB2,isect)
c*******************************************************************************
      use measurements, ONLY: xA1, yA1, xB1, yB1, nw, nh, n1, n2, n3, n4
     + , baseA, teeA, baseB, teeB
      real*8, dimension(100) :: xA2, yA2, xB2, yB2
      real*8, dimension(nw,*) :: topA, topB
      real*8, dimension(nh,*) :: sideA, sideB
c
      xA2 = xA1
      yA2 = yA1
      yA2(2:n1-1) = yA2(2:n1-1) + baseA(:,isect)
      xA2(n1+1:n2-1) = xA2(n1+1:n2-1) - teeA(nh:1:-1,isect)
      yA2(n2+1:n3-1) = yA2(n2+1:n3-1) + topA(nw:1:-1,isect)
      xA2(n3+1:n4-1) = xA2(n3+1:n4-1) - sideA(:,isect)
c
      xB2 = xB1
      yB2 = yB1
      yB2(2:n1-1) = yB2(2:n1-1) + baseB(:,isect)
      xB2(n1+1:n2-1) = xB2(n1+1:n2-1) + teeB(nh:1:-1,isect)
      yB2(n2+1:n3-1) = yB2(n2+1:n3-1) + topB(nw:1:-1,isect)
      xB2(n3+1:n4-1) = xB2(n3+1:n4-1) + sideB(:,isect)
c
c do corners
c
      xA2(1)=xA2(n4-1)
      xA2(n1)=xA2(n1+1)
      xA2(n2)=xA2(n2-1)
      xA2(n3)=xA2(n3+1)
      xA2(n4)=xA2(n4-1)
c
      yA2(1)=yA2(2)
      yA2(n1)=yA2(n1-1)
      yA2(n2)=yA2(n2+1)
      yA2(n3)=yA2(n3-1)
      yA2(n4)=yA2(2)
c
      xB2(1)=xB2(n4-1)
      xB2(n1)=xB2(n1+1)
      xB2(n2)=xB2(n2-1)
      xB2(n3)=xB2(n3+1)
      xB2(n4)=xB2(n4-1)
c
      yB2(1)=yB2(2)
      yB2(n1)=yB2(n1-1)
      yB2(n2)=yB2(n2+1)
      yB2(n3)=yB2(n3-1)
      yB2(n4)=yB2(2)
c
c Scale position along edges for cleaner plots (no overlapping corners)
c
      xA2(2:n1-1) = (/(xA2(1)+(xA2(n1)-xA2(1))*(i-.5)/nw,i=1,nw)/)
      yA2(n1+1:n2-1) = (/(yA2(n1)+(yA2(n2)-yA2(n1))*(i-.5)/nh,i=1,nh)/)
      xA2(n2+1:n3-1) = (/(xA2(n2)+(xA2(n3)-xA2(n2))*(i-.5)/nw,i=1,nw)/)
      yA2(n3+1:n4-1) = (/(yA2(n3)+(yA2(n4)-yA2(n3))*(i-.5)/nh,i=1,nh)/)
c
      xB2(2:n1-1) = (/(xB2(1)+(xB2(n1)-xB2(1))*(i-.5)/nw,i=1,nw)/)
      yB2(n1+1:n2-1) = (/(yB2(n1)+(yB2(n2)-yB2(n1))*(i-.5)/nh,i=1,nh)/)
      xB2(n2+1:n3-1) = (/(xB2(n2)+(xB2(n3)-xB2(n2))*(i-.5)/nw,i=1,nw)/)
      yB2(n3+1:n4-1) = (/(yB2(n3)+(yB2(n4)-yB2(n3))*(i-.5)/nh,i=1,nh)/)
c
      return
      end
c*******************************************************************************
      subroutine InitialSet
c*******************************************************************************
      use measurements
c calculate initial side clamp settings and top clamp targets
c
      write(*,"(' Initial Side Clamp Setting and Top Clamp Target')")
      write(*,"(' Clamp#  Side A  Side B  Top A   Top B')")
C      write(*,"('          -w1      w2      h1      h2 ')")
      write(*,"('         w-w1+x1 w+w2-x2   h1      h2 ')")
      do isect = 1,nsect
c current center targets
        dx = dxc(isect)
        dy = dyc(isect)
c Fit to curvature and torsion
        cinpl = Curv_InPl(isect)
        oop = Curv_OOP(isect)
        tor = Tors(isect)
        oopA = 1./(1./oop - dr)
        oopB = 1./(1./oop + dr)
        dh1= -0.0523757 + 0.1675188*abs(cinpl)
     +       + 0.7550988*abs(oopa) + 0.6649951*abs(tor)
        dh2= -0.0523757 + 0.1675188*abs(cinpl)
     +       + 0.7550988*abs(oopb) + 0.6649951*abs(tor)
        s1 = dh1*w
        s2 = dh2*w
c MCWF measurements
        x1 = -sum(teeA(1:nh,isect))/nh
        x2 =  sum(teeB(1:nh,isect))/nh
        y1 =  sum(baseA(1:nw,isect))/nw
        y2 =  sum(baseB(1:nw,isect))/nw
c Mike Z's area perserving fit based on initial estimates
        w1=-w/h*(y1+y2) - x2 - (s2+s1)/2/h + dx*2 + w/h*dy*2
        w2=+w/h*(y1+y2) - x1 + (s1+s2)/2/h + dx*2 - w/h*dy*2
        h1=-h/w*(x1+x2) - y2 - (s2-s1)/2/w + dy*2 + h/w*dx*2
        h2=+h/w*(x1+x2) - y1 + (s2-s1)/2/w + dy*2 - h/w*dx*2
        sideA5(:,isect)=-w1
        sideB5(:,isect)= w2
        topA5(:,isect) = h1
        topB5(:,isect) = h2
C      write(*,"(a4,1x,4f8.3)") name(isect), -w1, w2, h1, h2
      write(*,"(a4,1x,4f8.3)") name(isect),w-w1+x1, w+w2-x2, h1, h2
      enddo
c
      return
      end