c
c
c------------------------------------------------------------------------
c
c  main program for  f e r d o  allocates space for arrays
c
c------------------------------------------------------------------------
c
c***  program main ( input ,       output ,
c
c*** 1                             brespf ,        bwindf,
c
c*** 2               tape5=input , tape6=output,
c
c*** 3                             tape20=brespf , tape10=bwindf )
c
c
      dimension d(35000),id(35000)
c
      equivalence (d,id)
c
c------------------------------------------------------------------------
c
c  need length .ge. 21*nr + 5*nc + 13*nw + (nrc+nr)*nc
c
c                  + max(0,2*nloc-nr*nc,(nw-nrc-nr)*nc)
c
c   where:  nrc = max(nr,nc)
c
      length = 35000
c
      call setlo ( d, id, length )
c
      stop
c
      end 
c
c
c**********************************************************************
c
      subroutine setlo ( d, id, length )
c
c----------------------------------------------------------------------
c
c   reads dimensioning param's and sets length of variable arrays
c
c----------------------------------------------------------------------
c
      common /pltcon/ mi1,mi100,int100,int5,fl45,zero,izero,
c
     1                sevfiv,ten25,sixfiv,two75,fl90
c
      dimension jdent(20),d(length),id(length)
c
c-----------------------------------------------------------------------
c
      ierr = 0
c
c  set the value of the parameter tau used for biased estimation
c
      tau = 0.1
c
      read(5,501) jdent
      write(6,601) jdent
c
      read(5,502) nr,nc,nw,nloc
      write(6,602) nr,nc,nw,nloc
      nrc = max0(nr,nc)
c
      lblab = 1
      lb    = lblab + nr
      ls    = lb    + nr
      lut   = ls    + nr
      lllo  = lut   + 2*nr
      llup  = lllo  + nr
      lvlo  = llup  + nr
      lvup  = lvlo  + 2*nr
      lfort = lvup  + 2*nr
      lxplt = lfort + 2*nr
      lyplt = lxplt + 4*nr
      lw    = lyplt + 4*nr
      lq    = lw    + nc
      latmp = lq    + nc
      lelab = latmp + nc
      listr = lelab + nc
      lpest = listr + nc
      lplo  = lpest + nw
      lpup  = lplo  + nw
      lwv1  = lpup  + nw
      lwv2  = lwv1  + 2*nw
      lwv3  = lwv2  + 2*nw
      lwv4  = lwv3  + 2*nw
      lplab = lwv4  + 2*nw
      lwinw = lplab + nw
      la    = lwinw + nw
      lht   = la    + nrc*nc
      lxtra = lht   + nr*nc
      lnth  = lxtra + max0(0,2*nloc-nr*nc,nw*nc-2*nr*nc)
      lwdos = la
c  windos shares space with  a + ht
c  ht is used by subs ferdo and towcbn as working space
c
      write(6,603) length,lnth
      if(lnth.gt.length)  go to 99
c
   98 continue
c
      mi1 = -1
      mi100 = -100
      int100 = 100
      int5 = 5
      fl45 = 45.0
      zero = 0.0
      izero = 0
      sevfiv = 7.5
      ten25 = 10.25 
      sixfiv = 6.1
      two75 = 2.60
      fl90 = 90.0
c
      call qms2
c
c-------------------------------------------------------------------------
c
      call cntrlo (d(la),d(lb),d(ls),d(lw),d(lht),d(lq),d(lut),d(latmp),
c
     1  d(lplo),d(lpup),d(lblab),d(lplab),d(lwinw),d(lelab),id(listr),
c
     2  d(lxplt),d(lyplt),id(lllo),id(llup),d(lvlo),d(lvup),d(lfort),nr,
c
     3  nc,nw,nloc,tau,ierr,d(lwdos),
c
     4    d(lpest), d(lwv1), d(lwv2), d(lwv3), d(lwv4)  )
c
c------------------------------------------------------------------------
c
      call donepl
c
      return
c
   99 ierr = 1
      write(6,604)
      go to 98
c
  501 format(20a4)
  502 format(4i10)
  601 format(1h1,2x,20a4)
  602 format(1h0,' dimensioning parameters:',//,
     1 '     nr=',i5,3x,'no. of rows'/, 
     2 '     nc=',i5,3x,'no. of columns'/,
     3 '     nw=',i5,3x,'no. of windows'/,
     4 '   nloc=',i5,3x,'space needed for unbinned data'/)
  603 format(1h0,25(1h*),/i7,' locations allocated'/i7,
     1 ' locations required'/26(1h*)//) 
  604 format(1h0,' **** insufficient space - input options will be',
     1 ' edited, but problem will not be run'//)
c
c
      end 
c
c
c**********************************************************************
c
c    main control subroutine for  f e r d o
c
c----------------------------------------------------------------------
c
      subroutine cntrlo (a,b,s,w,ht,q,ut,atemp,plo,pup,blab,
c
     1                   plab,winw,elab,istar,xplwrk,yplwrk,llo,lup,
c
     2                   vlo,vup,forte,mr,mc,kr,nloc,tau,ierr,windos, 
c
     3                   pest, workv1, workv2, workv3, workv4 )
c
c----------------------------------------------------------------------
c
      real mu
c
      character*2 isym(3)
c
      dimension a(mr,mc),b(mr),s(mr),w(mc),ht(mr,mc),q(mc),ut(2*mr),
c
     1          atemp(mc),plo(kr),pup(kr)
c
      dimension  pest(kr), workv1(2*kr), workv2(2*kr), workv3(2*kr),
c
     1                     workv4(2*kr) 
c
      dimension windos(kr,mc),istar(mc) 
c
      dimension blab(mr),plab(kr),winw(kr),elab(mc)
c
      dimension ident(20),iaplt(7),iwplt(7)
c
      dimension  xplwrk(4*mr), yplwrk(4*mr)
c
      dimension llo(mr),lup(mr),vlo(2*mr),vup(2*mr),forte(2*mr)
c
      data nopt5,sk1,sk2/ 0, 0.0, 0.0/
c
      data isym / 'l ', 'u ', '= ' /
c
c----------------------------------------------------------------------
c
      nr = mr
      nc = mc
      nw = kr
c
      nbits = i1mach ( 11 )
c
      print *, ' '
      print *, ' nbits = ', nbits
      print *, ' '
c
c     nscr is logical unit no. for windows and a-matrix scratch
c
      nscr = 2
c
c
c  read window options card
   51 read(5,903) iwrd,iwwrt,strech,iwprt,nwplt
      write(6,921) iwrd,iwwrt,strech,iwprt,nwplt
      if(nwplt.eq.0)  go to 52
      read(5,901) (iwplt(i),i=1,nwplt)
      write(6,931) (iwplt(i),i=1,nwplt) 
c
c
c  read response matrix options card
   52 read(5,902) iard,iawrt,iaprt,naplt,nopt5,sk1,sk2
      write(6,922) iard,iawrt,iaprt,naplt,nopt5,sk1,sk2
      if(naplt.eq.0)  go to 53
      read(5,901) (iaplt(i),i=1,naplt)
      write(6,932) (iaplt(i),i=1,naplt) 
c
c
c  read measured data and output options card
c
   53 read(5,901) ibin,idiv,ngain,iovl,ibplt,ipplt,ippun,ihist
      write(6,923) ibin,idiv,ngain,iovl,ibplt,ipplt,ippun
c
c        stop if ierr \ 0  - insufficient space
      if(ierr.gt.0)  return
c
c  set value of mu for 99.95 pct. confidence in q-bounds (initial bounds)
c  This is not the value of mu used in the final bounds.
c
      twonm1 = 2*nr - 1
      mu = (1.0/sqrt(2.0))*(3.270 + sqrt(twonm1) )
      write(6,904) mu
c
c***  read plab, window functions, and elab
c
      if(iwrd.ge.0)  go to 21 
c    --------------------------------
      call matin ( plab, kr, nw, 1 )
c    --------------------------------
      call matin ( winw, kr, nw, 1 )
c    --------------------------------
      call matin ( elab, nc, nc, 1 )
c    --------------------------------
      do 20 i=1,nw
      winw(i) = winw(i) * strech
      call wread(w,nc,kr,i,elab,plab,winw)
      do 20 j=1,nc
   20 windos(i,j) = w(j)
      go to 23
c
   21 if(iwrd.gt.0)  go to 22 
      call matin(plab,kr,nw,1)
      call matin(windos,kr,nw,nc)
      call matin(elab,nc,nc,1)
      go to 23
c
   22 read( 10 ) plab
      read( 10 ) windos
      read( 10 ) elab
      if(iwrd.ne.iard)  rewind  10
c
c
c***  write window functions
c
   23 if(iwwrt.eq.0)  go to 24
      write(iwwrt) plab
      write(iwwrt) windos
      write(iwwrt) elab
      if(iwwrt.ne.iawrt)  rewind iwwrt
c
   24 do 26 i=1,nw
      do 25 j=1,nc
   25 w(j) = windos(i,j)
   26 write(nscr) w 
c
c  plot/print window functions as desired
      if(nwplt.eq.0)  go to 27
      call rsplot(windos,kr,mc,elab,nw,nc,nwplt,iwplt,1,xplwrk)
   27 if(iwprt.eq.0)  go to 28
      call vecpr ( plab, kr, 'plab', 1 )
      write(6,601)
      call wotac(windos,nw,nc,1,' row','col.','    ',6)
      call vecpr ( elab, nc, 'elab', 1 )
c
c***  read/write the matrix a 
c
   28 if(iard.eq.0)  call matin (a,mr,nr,nc)
c
      if ( iard .gt. 0 ) then 
                         read ( 20 ) a
                         rewind  20
      end if
c
c
   10 if(iawrt.eq.0) go to 11 
      write (iawrt) a
      rewind iawrt
c
   11 write(nscr) a 
      backspace nscr
c
c  plot/print chosen response functions 
      if(naplt.eq.0)  go to 12
      call rsplot(a,mr,mc,elab,nr,nc,naplt,iaplt,0,xplwrk)
   12 if(iaprt.eq.0)  go to 4 
      write(6,602)
      call wotac(a,nr,nc,1,' row','col.','    ',6)
c
c
c  program loops back to here for each new measured count vector (b)
c
    4 read (5,905,end=1040)   ident
      write (6,906)  ident
c
c  If ibin.eq.0 (i.e., if data already binned), the blab, b, and s
c    are read by matin in bsread.  otherwise, towcbn, et al., will
c    read channel data and calc blab, b, and s.  towcbn needs ht
c    for working space.
c
c    ----------------------------------------------------------
      call bsread ( b,s,mr,nr,nc,blab,ibin,idiv,ngain,llo,lup,
c
     1              vlo,vup,forte,iovl,ht,2*nloc )
c    ----------------------------------------------------------
c
c  Read numbers of first and last window responses to include in output plots
c
      read ( 5, * ) n1, n2
c
      nwrpl = n2 - n1 + 1
c
      print *, ' '
      print *, ' no. of first window response to plot = ', n1
      print *, ' no. of last  window response to plot = ', n2
c
      if(naplt.eq.0)  go to 14
c
      call arplot ( a,mr,mc,blab,nr,nc,naplt,iaplt,0,xplwrk )
c
      naplt = 0
c
   14 continue
      if(ibplt.eq.0.or.ibplt.eq.2)  go to 308
      write(6,910) ident
c
      call aplot ( b,mr,1,2,nr,1,2,blab,dummy )
c
  308 if(ibplt.gt.1)
c                  -------------------------------------------------
     1              call plotb ( b,s,mr,nr,blab,xplwrk,yplwrk,4*mr,
c
     2                           forte, vlo, vup, ut )
c                  -------------------------------------------------
c
      write(6,906) ident
c
c---------------------------------------------------------------------
c
      call ferdo ( a,b,s,w,ht,q,ut,atemp,istar,plab,mr,mc,
c
     1             nr,nc,nw,nbits,tau,kr,mu,plo,pup,nscr,
c
     2             pest )
c
c---------------------------------------------------------------------
c
      if(ippun.eq.0)  go to 6 
      write(7,905) ident
      write(7,901) nw
      write(7,909) (plab(i),plo(i),pup(i),i=1,nw) 
c
    6 if(ipplt.eq.0.or.ipplt.eq.2)  go to 5
      write (6,907) ident
c
      call pplot3 (plo,nw,6,nw,0.004,2,4hener,plab,isym,6)
c
      write (6,908) ident
c
      call aplot(plo(6),kr,1,2,nw-5,1,2,plab(6),6)
c
    5 if(ipplt.lt.2)  go to 4 
c
      if (ihist.eq.1) then
c          ------------------------------------------------------ 
            call phplot ( plo, pup, n1, nw, nwrpl, kr, plab,
c
     1                    xplwrk, yplwrk, 4*mr,
c
     1                    pest, workv1, workv2, workv3, workv4,
c
     2                    'FERDO Estimates $' ) 
c          ------------------------------------------------------
            else
c            -----------------------------------------------------
              call phplot ( plo, pup, n1, nw, nwrpl, kr, plab,
c
     1                      xplwrk, yplwrk, 2*mr,
c
     1                      pest, workv1, workv2, workv3, workv4,
c
     2                      'FERDO Estimates $' ) 
c            -----------------------------------------------------
      end if
c
c
c    go back to read in another count vector (b)
c
      go to 4
c
c
 1040 return
c
c
  601 format(1h0,//' window matrix'/)
  602 format(1h0,//' response matrix'/) 
  901 format(8i10)
  902 format(5i10,2e10.0)
  903 format (2i10,f10.3,2i10)
  904 format(1h0/1h0,'*****   mu  =',1pe11.3,'   *****'/)
  905 format(20a4)
  906 format(1h1,20a4)
  907 format(1h1,'linear plot of plo and pup',6x,20a4)
  908 format(1h1,'pseudo-log plot of plo and pup',6x,20a4)
  909 format(6e12.5)
  910 format(1h1,'pseudo-log plot of b and s',6x,20a4)
  921 format(1h0,' window options:'//,
     1 '   iwrd=',i5,3x,'-1 - calc guassn windows from plab + winw'/, 
     *              16x,' 0 - read 2d window matrix from cards'/,
     *              16x,' n - read plab, windos, + elab from unit n'/,
     2 '  iwwrt=',i5,3x,'0/n - no effect/write plab, windos, and elab',
     * ' on unit n'/,
     3 ' strech=',f5.2,3x,'mult. factor for winw (if iwrd=-1)'/,
     4 '  iwprt=',i5,3x,'0/1 - no effect/print 2d window matrix'/,
     5 '  nwplt=',i5,3x,'0/n - no. of window fcts to plot (nwplt)8)') 
  922 format(1h0,/'  response matrix options:'//, 
     1 '   iard=',i5,3x,' 0 - read response matrix from cards'/,
     *              16x,' n - read response matrix from unit n'/,
     2 '  iawrt=',i5,3x,'0/n - no effect/write response on unit n'/,
     3 '  iaprt=',i5,3x,'0/1 - no effect/print response matrix'/,
     4 '  naplt=',i5,3x,'0/n - no. of resp fcts to plot (naplt)8)'/,
     5 '  nopt5=',i5,3x,'not used in ferdo'/,
     6 '    sk1=',f5.2,3x,'not used in ferdo'/,
     7 '    sk2=',f5.2,3x,'not used in ferdo')
  923 format(1h0,/'  measured data and output options:'//,
     1 '   ibin=',i5,3x,'0/1 - measured data has/hasnot been binned'/,
     1 '   idiv=',i5,3x,'0/1 - no effect/divide binned data by ',
     & 'bin width (if ibin = 1)',/,
     2 '  ngain=',i5,3x,'no. of overlapping gains in data (ngain)3)'/,
     3 '   iovl=',i5,3x,'first bin of 4-bin overlap (ngain=2)'/,
     4 '  ibplt=',i5,3x,'0 - no effect'/,
     *              16x,'1 - printer plot of binned spectrum'/,
     *              16x,'2 - calcomp plot of binned spectrum'/,
     *              16x,'3 - both printer and calcomp plots'/,
     5 '  ipplt=',i5,3x,'same as ibplt except for plo-pup solution'/, 
     6 '  ippun=',i5,3x,'0/1 - no effect/punch plo-pup solution'/)
  931 format(1h0,' window functions to plot:',7i6)
  932 format(1h0,' response functions to plot:',7i6)
      end 
