      program cal_search17_2
***********************************************************************
*     search calibration files                                        *
*                                    October 1, 2008  H. Koshiishi    *
***********************************************************************
c
      integer naxes(2),nantst(84),icimg(1000),iimg(1000)
      character*80 filehd,infpar,outf
      character*12 dateframe
c
ccc   write(6,'(/'' input directory for output files       : ''$)')
      read(5,'(a)') outf
ccc   write(6,'('' input rawdata filename header          : ''$)')
      read(5,'(a)') filehd
ccc   write(6,'('' input parameter filename               : ''$)')
      read(5,'(a)') infpar
c
ccc   write(6,'('' input first frame number               : ''$)')
      read(5,*) ncfrst
      if (ncfrst.ne.0) then
ccc    write(6,'('' input last frame number                : ''$)')
       read(5,*) ncfrend
ccc    write(6,'('' input frame interval                   : ''$)')
       read(5,*) ictvl
       ncimg=0
       do 150 i=ncfrst,ncfrend,ictvl
        ncimg=ncimg+1
        icimg(ncimg)=i
  150  continue
      else
ccc    write(6,'('' input number of frames                 : ''$)')
       read(5,*) ncimg
ccc    write(6,'('' input frame numbers                    : ''$)')
       read(5,*) (icimg(i),i=1,ncimg)
      end if
ccc   write(6,'('' input integration frames'')')
ccc   write(6,'(''      for calibration                   : ''$)')
      read(5,*) iccal
ccc   write(6,'(''      for image processing              : ''$)')
      read(5,*) ictgr
c
ccc   write(6,'('' quiet : 1 / burst : 2               : ''$)')
      read(5,*) iqb
c
ccc   write(6,'('' input pixel-size of the output files   : ''$)')
      read(5,*) npix
c
ccc   write(6,'('' input offset of the image center'')')
ccc   write(6,'('' unit(1:arcsec/2:pixel) & x-/y-offsets  : ''$)')
      read(5,*) nuoff,xoffa,yoffa
      if (nuoff .eq. 1) then
       xoffp=xoffa/4.91106
       yoffp=yoffa/4.91106
      else
       xoffp=xoffa
       yoffp=yoffa
       xoffa=xoffp*4.91106
       yoffa=yoffp*4.91106
      endif
c
ccc   write(6,'('' input first frame number               : ''$)')
      read(5,*) nfrst
      if (nfrst.ne.0) then
ccc    write(6,'('' input last frame number                : ''$)')
       read(5,*) nfrend
ccc    write(6,'('' input frame interval                   : ''$)')
       read(5,*) itvl
       nimg=0
       do 160 i=nfrst,nfrend,itvl
        nimg=nimg+1
        iimg(nimg)=i
  160  continue
      else
ccc    write(6,'('' input number of frames                 : ''$)')
       read(5,*) nimg
ccc    write(6,'('' input frame numbers                    : ''$)')
       read(5,*) (iimg(i),i=1,nimg)
      end if
ccc   write(6,'('' input integration frames'')')
ccc   write(6,'(''      for calibration                   : ''$)')
      read(5,*) ical
ccc   write(6,'(''      for image processing              : ''$)')
      read(5,*) itgr
c
ccc   write(6,'('' input number limit of criterion(r+l,r-l): ''$)')
      read(5,*) nmaxa,nmaxs
ccc   write(6,'('' input lower limit of criterion(r+l,r-l): ''$)')
      read(5,*) cradd,crsub
ccc   write(6,'('' input factor of criterion              : ''$)')
      read(5,*) crfac
c
***** default parameters *****
c
      loutf=index(outf,' ')-1
c
      call cft512(mapra,maprb,0)
c
      iunit=11
c
      iounit=2
      nbit=32
      naxis=2
      naxes(1)=336
      naxes(2)=1
c
      do 100 iifr=1,ncimg
c
      iframe=icimg(iifr)
c
***** calibration *****
c
      call rddcal17_2(iunit,filehd,iframe,iccal,nantst,
     -     phnsr,phewr,phnsl,phewl)
c
      lfilehd=index(filehd,' ')-1
      write(dateframe,'(a6,''_'',i5)')filehd(lfilehd-5:lfilehd),iframe
c
      phr=sqrt((phnsr*phnsr+phewr*phewr)/2.)
      phl=sqrt((phnsl*phnsl+phewl*phewl)/2.)
c
      write(6,'(a12,x,a12)') 'date_frame =',dateframe
      write(6,'(a24,x,f8.2,x,f8.2,x,f8.2)')
     - 'phns , phew , ph (rcp) =',phnsr,phewr,phr
      write(6,'(a24,x,f8.2,x,f8.2,x,f8.2)')
     - 'phns , phew , ph (lcp) =',phnsl,phewl,phl
      write(6,*) ' '
c
  100 continue
c
      end
c
c
***********************************************************************
*     rddcal17_2 - read raw heliograph data and calibrate             *
*       iunit  : input,integer fortran file unit number               *
*                  for raw data file                                  *
*       filehd : input,character*80 raw data filename header          *
*       iframe : input,integer  frame number for calibration          *
*       itgr   : input,integer  number of frames for integration      *
*       nantst : input,integer(84) antenna status                     *
*       sums_phansr : rcp/ph/ns-cal error                             *
*       sums_phaewr : rcp/ph/ew-cal error                             *
*       sums_phansl : lcp/ph/ns-cal error                             *
*       sums_phaewl : lcp/ph/ew-cal error                             *
*                                    October 2, 2008  H. Koshiishi    *
***********************************************************************
c
      subroutine rddcal17_2(iunit,filehd,iframe,itgr,nantst,
     -     sums_phansr,sums_phaewr,sums_phansl,sums_phaewl)
c
      character*80 filehd
      integer nantst(84)
      real corrc(84,85),corrs(84,85),corlc(84,85),corls(84,85)
c
      integer iposns(28),iposew(57),icv(9)
      real*8 pha(200),phansr(2000),phansl(2000),
     -       phaewr(2000),phaewl(2000)
c
      integer nvco(2000),nac1(2000),nac2(2000)
      real*8 amm(2000,200),amlo(2000),amo(2000)
c
      data iposns/0,1,2,3,4,5,6,7,8,9,10,11,
     -          12,14,16,20,24,28,
     -          32,40,48,56,
     -          64,80,96,112,128,144/
      data iposew/-160,-144,-128,-112,-96,-80,-64,
     -          -56,-48,-40,-32,
     -          -28,-24,-20,-16,-14,-12,
     -          -11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,
     -          0,1,2,3,4,5,6,7,8,9,10,11,
     -          12,14,16,20,24,28,
     -          32,40,48,56,
     -          64,80,96,112,128,144,160/
      data icv/1,2,4,6,8,12,16,24,32/
c
      nantns=28
      nantew=57
      ncv=9
c
      call rddat17(iunit,filehd,iframe-(itgr-1)/2,
     -     iframe+itgr/2,nantst,
     -     ndjst,msjst,jsts,jste,natt,nfalt,nfreq,nstat,
     -     corrc,corrs,corlc,corls,nframe)
c
      inic=0
c
      call ammat(nantns,iposns,ncv,icv,nrow,nvco,nac1,nac2,amm)
      call amobs(corrc,corrs,1,nrow,nac1,nac2,amo,amlo)
      call phcal_2(nantns,iposns,ncv,icv,1,corrc,corrs,inic,pha,phansr)
      sums_phansr=0.
      num=0.
      do 100 il=1,nrow
         if (amo(il).ne.0.) then
            phansr(il)=phansr(il)/amo(il)*1000.
            sums_phansr=sums_phansr+phansr(il)*phansr(il)
            num=num+1.
         endif
  100 continue
      sums_phansr=sqrt(sums_phansr/num)
c
      call ammat(nantew,iposew,ncv,icv,nrow,nvco,nac1,nac2,amm)
      call amobs(corrc,corrs,2,nrow,nac1,nac2,amo,amlo)
      call phcal_2(nantew,iposew,ncv,icv,2,corrc,corrs,inic,pha,phaewr)
      sums_phaewr=0.
      num=0.
      do 200 il=1,nrow
         if (amo(il).ne.0.) then
            phaewr(il)=phaewr(il)/amo(il)*1000.
            sums_phaewr=sums_phaewr+phaewr(il)*phaewr(il)
            num=num+1.
         endif
  200 continue
      sums_phaewr=sqrt(sums_phaewr/num)
c
      call ammat(nantns,iposns,ncv,icv,nrow,nvco,nac1,nac2,amm)
      call amobs(corlc,corls,1,nrow,nac1,nac2,amo,amlo)
      call phcal_2(nantns,iposns,ncv,icv,1,corlc,corls,inic,pha,phansl)
      sums_phansl=0.
      num=0.
      do 300 il=1,nrow
         if (amo(il).ne.0.) then
            phansl(il)=phansl(il)/amo(il)*1000.
            sums_phansl=sums_phansl+phansl(il)*phansl(il)
            num=num+1.
         endif
  300 continue
      sums_phansl=sqrt(sums_phansl/num)
c
      call ammat(nantew,iposew,ncv,icv,nrow,nvco,nac1,nac2,amm)
      call amobs(corlc,corls,2,nrow,nac1,nac2,amo,amlo)
      call phcal_2(nantew,iposew,ncv,icv,2,corlc,corls,inic,pha,phaewl)
      sums_phaewl=0.
      num=0.
      do 400 il=1,nrow
         if (amo(il).ne.0.) then
            phaewl(il)=phaewl(il)/amo(il)*1000.
            sums_phaewl=sums_phaewl+phaewl(il)*phaewl(il)
            num=num+1.
         endif
  400 continue
      sums_phaewl=sqrt(sums_phaewl/num)
c
      return
      end
c
c
***********************************************************************
*     phcal_2 : phase calibration subroutine                          *
*       assumption : phase s1 = 0, a0 = 0                             *
*       using a part of correlations                                  *
*       weighted errors in least square method                        *
*       double precision                                              *
*                                                                     *
*     nant : input,integer  number of antennas                        *
*     ipos : input,integer  antenna positions                         *
*       dimension of ipos : number of antennas (nant)                 *
*     ncv  : input,integer  number of correlation vector              *
*     icv  : input,fourier  component vector                          *
*       dimension of icv : number of different vector (ncv)           *
*     nsew : input,integer  (1)ns or (2)ew                            *
*     corc : input,real     cos of correlation                        *
*     cors : input,real     sin of correlation                        *
*     inic : input,integer  condition initial value                   *
*     pha  : input,real     initial values if inic ne 0               *
*     pha  : output,real    antenna phases                            *
*     phc  : output,real    phase errors                              *
*                                                                     *
*     phm : phase matrix                                              *
*     amo : observed correlation amplitudes                           *
*     pho : observed correlation phases                               *
*     phc : calculated errors of correlation phases                   *
*     phx : calculated antenna phases and fourier phases              *
*     vw : work area                                                  *
*                                    October 15,2008  H. Koshiishi    *
***********************************************************************
c
      subroutine phcal_2(nant,ipos,ncv,icv,nsew,corc,cors,inic,
     -     pha,phc)
c
      implicit real*8 (a-h,o-z)
      integer ipos(1),icv(1),erst(100),nvco(2000),nac1(2000),
     -     nac2(2000)
      real*8 phm(2000,200),pho(2000),amo(2000),phc(2000),
     -     phx(200),vw(2200,200),ang(50)
      real corc(84,85),cors(84,85),pha(1),rerr
c
      external pherr
c
      common /phmatr/nrow,ncol,phm,amo,pho,nac1,nac2
c
***** matrix for calibration *****
c
      ncol=nant-1+ncv-1
      call phmat(nant,ipos,ncv,icv,nrow,nrow1,nrow2,nvco,nac1,nac2,phm)
c
***** observed correlation phases *****
c
      call phobs(corc,cors,nsew,nrow,nac1,nac2,amo,pho)
c
***** initial estimation *****
c
      do 400 i=1,nant
  400 if (ipos(i).eq.0) ncent=i
c
      if (inic.eq.0) then
c
      do 420 i=1,nant
  420 erst(i)=0
      do 440 ic=1,ncol
  440 phx(ic)=0.
c
*     initial estimation for east antennas
c
      if (nrow2.gt.nrow1) then
      iv=0
      do 610 ir=nrow1+1,nrow2
       do 620 i=nant-1,nant-ncent+1,-1
  620  if (phm(ir,i).eq.-1.) go to 630
  630  iant2=i
       do 640 i=iant2,nant-ncent+1,-1
  640  if (phm(ir,i).eq.1.) go to 650
       i=0
  650  iant1=i
       if (nvco(ir)-1.ne.iv) nsamp=0
       iv=nvco(ir)-1
c
       if (erst(iant2).eq.0) then
        if (iv.eq.0) then
         if (iant1.eq.0) then
          phx(iant2)=-pho(ir)
         else
          phx(iant2)=phx(iant1)-pho(ir)
         end if
        else
         if (iant1.eq.0) then
          phx(iant2)=phx(nant-1+iv)-pho(ir)
         else
          phx(iant2)=phx(iant1)+phx(nant-1+iv)-pho(ir)
         end if
        end if
        phx(iant2)=dmod(dmod(phx(iant2)+180.,360.d0)+360.,360.d0)-180.
        erst(iant2)=1
       else
        nsamp=nsamp+1
        if (iant1.eq.0) then
         ang(nsamp)=pho(ir)+phx(iant2)
        else
         ang(nsamp)=pho(ir)-phx(iant1)+phx(iant2)
        end if
c        write(6,'('' '',i5,5f10.2)') nsamp,(ang(isamp),isamp=1,nsamp)
        call dangmn(nsamp,ang,phx(nant-1+iv))
       end if
  610 continue
      do 670 i=nant,ncol
  670 phx(i)=-phx(i)
c
      end if
c
*     initial estimation for north or west antennas
c
      if (nrow1.ne.0) then
c
      iv=0
      do 500 ic=nant,ncol
  500 phx(ic)=0.
      do 510 ir=1,nrow1
       do 520 i=nant-ncent,1,-1
  520  if (phm(ir,i).eq.-1.) go to 530
  530  iant2=i
       do 540 i=iant2,1,-1
  540  if (phm(ir,i).eq.1.) go to 550
       i=0
  550  iant1=i
       if (nvco(ir)-1.ne.iv) nsamp=0
       iv=nvco(ir)-1
c
       if (erst(iant2).eq.0) then
        if (iv.eq.0) then
         if (iant1.eq.0) then
          phx(iant2)=-pho(ir)
         else
          phx(iant2)=phx(iant1)-pho(ir)
         end if
        else
         if (iant1.eq.0) then
          phx(iant2)=phx(nant-1+iv)-pho(ir)
         else
          phx(iant2)=phx(iant1)+phx(nant-1+iv)-pho(ir)
         end if
        end if
        phx(iant2)=dmod(dmod(phx(iant2)+180.,360.d0)+360.,360.d0)-180.
        erst(iant2)=1
       else
        nsamp=nsamp+1
        if (iant1.eq.0) then
         ang(nsamp)=pho(ir)+phx(iant2)
        else
         ang(nsamp)=pho(ir)-phx(iant1)+phx(iant2)
        end if
        call dangmn(nsamp,ang,phx(nant-1+iv))
       end if
  510 continue
c
      end if
c
      else
      do 580 i=1,ncol
  580 phx(i)=pha(i)
      end if
c
***** least square calculation *****
c
      nmax=nrow*ncol*100
      epsr=1.d-3
      call dnolf1(phx,ncol,pherr,nrow,epsr,nmax,
     -           phc,sums,vw,2200,icon)
      do 700 i=1,ncol
  700 pha(i)=dmod(dmod(phx(i)+180.,360.d0)+360.,360.d0)-180.
      rerr=sums
c
      if (icon.ne.0) write(6,'(''  phcal : icon of dnolf1 ='',i6)')
     -     icon
c
      return
      end
