subroutine mtcarlo(modvis2,modphas,mcinit,mcrun)
c
c     Analyze calibrator data and create artificial data
c
c     C. A. Hummel      10 November 1992
c
      save pdata,pdates,ptimmin,ptimstep,gonly,test,
     &     dopower,doboot,doshuffle,
     &     caldata,visdata,ndates,ncals,idum
c
      include 'constnts.inc'
      include 'data.inc'
      include 'model.inc'
c
      real*8 modvis2(mxrec,mxfilt),modphas(mxrec,mxfilt)
      real*8 visdata(mxrec,mxfilt)
      real*4 y(mxpl),t(mxpl),fpower(mxpl,mxfilt),fphase(mxpl,mxfilt)
      real*4 caldata(mxdates,mxstars,mxfilt,mxpl)
      real*4 timmax,timmin,timstep
      integer*4 firstid,lastid
      integer*4 boot1(mxrec),boot2(mxrec),ncals(mxdates)
      logical modplot,query,mcinit,mcrun,gonly,test
      logical doboot,dopower,doshuffle,firstsim,ask
      character header*80,device*4,plotfile*15,command*1

c     These are for the linear prediction programs
      integer*4 npoles
      parameter (npoles=5)
      real*4 my(mxpl),wk1(mxpl),wk2(mxpl),wkm(npoles),d(npoles)
      real*4 future(mxpl)

c     These are for the Fourier-transform
      real*4 data(2*mxpl),pdata(mxpl,mxfilt,mxdates),phi(mxpl)
      integer*4 pdates(mxdates)
      real*4 ptimmin(mxdates),ptimstep(mxdates)
c
      data itest/0/

      if(.not.mcinit)goto 30
      idum=-1
      write(outc,'(a,$)')'Restart? Enter random numbers to skip: '
      read(inc,*)ir
      do i=1,ir       
       x=ran1(idum)  ! Initialize random number generator
      enddo
c
c     Initialize
c
      if(query('Include calibration errors ? (y/n) '))then
       gonly=.false.
       dopower=.false.
       doboot=.false.
       doshuffle=.false.
       firstsim=.true.
       write(outc,*)
       write(outc,*)'Please select Simulation method: '
       write(outc,*)'-1- Power spectra'
       write(outc,*)'-2- Bootstrapping'
       write(outc,*)'-3- Calibrator shuffling'
       write(outc,'(a,$)')' >'
       read(inc,2)device
       if(device.eq.'1')then
        write(outdat,'(a)')'! MT CARLO: Power spectra method selected'
        dopower=.true.
       endif
       if(device.eq.'2')then
        write(outdat,'(a)')'! MT CARLO: Bootstrapping method selected'
        doboot=.true.
        do id=1,ndata       ! save observed visibility errors
         do ig=1,nfilt
          visdata(id,ig)=viserr(id,ig)
         enddo
        enddo
        mcinit=.false.
        return
       endif
       if(device.eq.'3')then
        write(outdat,'(a)')'! MT CARLO: Calibrator shuffling selected'
        doshuffle=.true.
       endif
      else
       write(outdat,'(a)')'! MT CARLO: Gaussian noise only  selected'
       gonly=.true.
       mcinit=.false.
       return
      endif
c
c     Initialize plotting                
c
      write(outc,*)'Screen (s) or hardcopy (h) or exit (n) ?'
      read(inc,2)device
    2 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotmtc.lpm'//'/lpm'
       iy=4
       csize=2.5
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
       iy=3
       csize=2.5
      else
       return
      endif
      ask=.true.
      if(query('Fast graphics ? (y/n) '))ask=.false.
c
      call pgbegin(0,plotfile,3,iy)
      call pgask(ask)
      call pgsch(csize)
      firstid=1
      lastid=1
      idate=1
      do ig=1,nfilt
       do i=1,mxpl*2,2
        pdata(i/2+1,ig,idate)=0.
       enddo
      enddo
c
c     Go through the whole data file and ...
      DO ID=1,NDATA
c
c      ... do analysis, when date or baseline changes
       if(date(id+1).ne.date(id).or.baseline(id+1).ne.baseline(id)
     &                          .or.id+1.gt.ndata)then
        timmin=real(hours(firstid))
        timmax=real(hours(lastid))
        timstep=(timmax-timmin)/real(mxpl-1)
        pdates(idate)=date(id)
        ptimmin(idate)=timmin
        ptimstep(idate)=timstep
        ndates=idate
c
c       Do one calibrator after the other
        ic=1
        ncal=0
        do while(calstars(ic).ne.0)
c
c ...... Find out when valid data (3 chs) for current star begins and ends
         ttot=24.
         do ig=1,nfilt
          tmin=0
          tmax=0
          do n=1,lastid-firstid+1
           if(stars(n+firstid-1).eq.calstars(ic))then
            if(viserr(n+firstid-1,ig).ne.2.D0)then
             tmin=real(hours(n+firstid-1))
             goto 10
            endif
           endif
          enddo
   10     do n=lastid-firstid+1,1,-1
           if(stars(n+firstid-1).eq.calstars(ic))then
            if(viserr(n+firstid-1,ig).ne.2.D0)then
             tmax=real(hours(n+firstid-1))
             goto 20
            endif
           endif
          enddo
   20     ttot=min(ttot,tmax-tmin)
         enddo
         if(ttot.gt.2)then   ! do only stars with more than 3h data
          ncal=ncal+1
c
c         loop over nfilt plots of residuals for given calibrator star
          do ig=1,nfilt
           ymax=0.
           ymin=0.
c
c ........ Do a running average and plot residuals at the same time
           do k=1,mxpl
            t(k)=timmin+timstep*(k-1)
            y(k)=0.
            swgn=0.
            if(k.eq.mxpl)then   ! Prepare for plot of y4,t4 on last run 
             call pgsci(5) 
             call pgenv(timmin,timmax,ymin,ymax,0,1) 
             write(header,3) 
     &        calstars(ic),date(id),baseline(id),int(lambda0(id,ig)) 
    3        format('FK',i4,1x,i6,1x,' B.',i2,1x,i3,' nm') 
             csize=csize+1 
             call pgsch(csize) 
             call pglabel('','',header) 
             csize=csize-1 
             call pgsch(csize) 
             if(ig.eq.1)call pgsci(2) 
             if(ig.eq.2)call pgsci(3) 
             if(ig.eq.3)call pgsci(11) 
            endif
            do n=1,lastid-firstid+1
             if(stars(n+firstid-1).eq.calstars(ic))then
              y4=real(obsvis2(n+firstid-1,ig)
     &               /modvis2(n+firstid-1,ig))-1.
              s4=real(viserr(n+firstid-1,ig))
              t4=real(hours(n+firstid-1))
              if(s4.ne.2.0)then
               ymax=max(ymax,y4)
               ymin=min(ymin,y4)
               if(abs(t4-t(k)).lt.3.0)then    
                wg=exp(-(t4-t(k))**2)         ! Gaussian weighting
                wn=1./(s4**2)                 ! Noise
                y(k)=y(k)+y4*(wg*wn)
                swgn=swgn+(wg*wn)
               endif
               if(k.eq.mxpl)then              ! Plot y4,t4 on last run
                call pgpoint(1,t4,y4,1)
                obsvis2lo=y4-s4
                obsvis2hi=y4+s4
                call pgerry(1,t4,obsvis2lo,obsvis2hi,1.0)
               endif
              endif
             endif
            enddo
            if(swgn.ne.0)then
             y(k)=y(k)/swgn
            else
             y(k)=0.
            endif
           enddo
c
c ........ Copy relevant running average points into array my
           nn=0
           do k=1,mxpl
            if(t(k).ge.tmin.and.t(k).le.tmax)then
             nn=nn+1
             my(nn)=y(k)
             noff=k-nn
            endif
           enddo
           nfut=(timmax-tmax)/timstep+1
c
c ........ Now extrapolate running average forward to timmax
           call memcof(my,nn,npoles,dum,d,wk1,wk2,wkm)
           call predic(my,nn,d,npoles,future,nfut)
           do n=1,nfut
            y(n+nn+noff)=future(n)
           enddo
c
c ........ Copy relevant running average points in reverse order in to my
           nn=0 
           do k=mxpl,1,-1 
            if(t(k).ge.tmin.and.t(k).le.tmax)then 
             nn=nn+1 
             my(nn)=y(k)
             noff=k+nn 
            endif 
           enddo 
           nfut=(tmin-timmin)/timstep+1 
c 
c ........ Now extrapolate running average backward to timmin 
           call memcof(my,nn,npoles,dum,d,wk1,wk2,wkm) 
           call predic(my,nn,d,npoles,future,nfut) 
           n=1
           do while(noff-nn-n.gt.0)
            y(noff-nn-n)=future(n)
            n=n+1
           enddo
c
c ........ Plot new array y
           call pgsci(6)
           call pgline(mxpl,t,y)
c
c......... Fill in array caldata
           do i=1,mxpl
            caldata(idate,ncal,ig,i)=y(i)
           enddo
c
c ........ Do FFT and power spectrum
           do i=1,mxpl*2,2
            data(i)=y(i/2+1)
            data(i+1)=0.
           enddo
           isign=1
           call four1(data,mxpl,isign)
           do i=1,mxpl*2,2
            pdata(i/2+1,ig,idate)=pdata(i/2+1,ig,idate)
     &                           +(data(i)**2+data(i+1)**2)
            fpower(i/2+1,ig)=data(i)**2+data(i+1)**2
            fphase(i/2+1,ig)=atan2(data(i),data(i+1))
           enddo
c
          enddo                ! Next filter
c
          if(dopower)then
           ymin=0.             ! Plot fourier phases
           ymax=0.
           do ig=1,nfilt
            do i=1,mxpl
             wk1(i)=real(i)
             my(i)=fphase(i,ig)/real(degrad)
             ymin=min(ymin,my(i))
             ymax=max(ymax,my(i))
            enddo
            call pgsci(5)
            call pgenv(0.,wk1(mxpl),ymin,ymax,0,0)
            call pglabel('Bin','Phase [degrees]','Fouriertransform')
            call pgsci(6)
            call pgline(mxpl,wk1,my)
           enddo
           ymin=0.             ! Plot powerspectrum amplitudes
           ymax=0.
           do ig=1,nfilt
            do i=1,mxpl
             wk1(i)=real(i)
             my(i)=log(fpower(i,ig))
             ymin=min(ymin,my(i))
             ymax=max(ymax,my(i))
            enddo
            call pgsci(5)
            call pgenv(0.,wk1(mxpl),ymin,ymax,0,0)
            call pglabel('Bin','log(Power)','Fouriertransform')
            call pgsci(6)
            call pgline(mxpl,wk1,my)
           enddo
          endif

         endif                 ! Not enough data for linear prediction
c
         ic=ic+1
        enddo                  ! Next calibrator
c
c ..... Take the mean of the power spectra
        if(ncal.eq.0)then
         idate=idate-1
        else
         do ig=1,nfilt
          do i=1,mxpl*2,2
           pdata(i/2+1,ig,idate)=sqrt(pdata(i/2+1,ig,idate)/ncal)
          enddo
         enddo
        endif
c
c ..... Now prepare for the next date
        ncals(idate)=ncal
        firstid=lastid+1
        lastid=firstid-1
        idate=idate+1
        do ig=1,nfilt
         do i=1,mxpl*2,2
          pdata(i/2+1,ig,idate)=0.
         enddo
        enddo
c
c .... Continue here, if date did not change
       endif
       lastid=lastid+1
      ENDDO                    ! End of ID-loop
      call pgask(.true.)
      call pgend
      mcinit=.false.
      return
c-----------------------------------------------------------------------
c     Now that we have initialized Monte Carlo power spectra,
c     we jump into the following code to simulate data.
c
   30 continue
      if(.not.mcrun.and.itest.eq.0)then
       if(query('View simulated data ? (y/n) '))then
        test=.true.
        call pgbegin(0,'/ega',3,3)
        device='s'
        modplot=.true.
       else
        test=.false.
       endif
      endif
      if(doshuffle.and.firstsim)then
       do id=1,ndata       ! save observed visibility amplitudes
        do ig=1,nfilt
         visdata(id,ig)=obsvis2(id,ig)
        enddo
       enddo
      endif
      firstsim=.false.
c.......................................
c     For gaussian noise only, do the following
      IF(gonly)THEN
       do ig=1,nfilt
        do id=1,ndata
         if(viserr(id,ig).ne.2.D0)
     &   obsvis2(id,ig)=modvis2(id,ig)
     &                 +dble(gasdev(idum))*viserr(id,ig)
        enddo
       enddo
      ELSE 
c.......................................
       if(doboot)then
c
        do id=1,ndata    
         boot2(id)=id
        enddo
        do id=1,ndata     
         x=ran1(idum)
         iid=int(x*ndata+1.)
         boot1(id)=boot2(iid)
         x=ran1(idum)
         iiid=int(x*ndata+1.) 
         boot2(iid)=iiid 
        enddo
        do id=1,ndata
         iweight=0
         do iid=1,ndata
          if(boot1(iid).eq.id)iweight=iweight+1
         enddo
         do ig=1,nfilt
          if(iweight.ne.0)then
           viserr(id,ig)=visdata(id,ig)/sqrt(dble(iweight))
           if(visdata(id,ig).eq.2.0D0)viserr(id,ig)=2.0D0 
          else
           viserr(id,ig)=2.0D0
          endif
         enddo
        enddo
c
       else
c
        firstid=1
        lastid=1
        DO ID=1,NDATA
        if(date(id+1).ne.date(id).or.baseline(id+1).ne.baseline(id)
     &                           .or.id+1.gt.ndata)then
         timmin=real(hours(firstid))
         timmax=real(hours(lastid))
         timstep=(timmax-timmin)/(mxpl-1)

         if(dopower)then
          idate=0
          do i=1,mxdates
           if(date(id).eq.pdates(i))idate=i
          enddo
          if(idate.eq.0)then     
           write(outc,*)
     &     '***ERROR: Did not process any calibrator for date ',date(id)
           goto 40
          endif
          do i=1,mxpl/2
           phi(i)=ran1(idum)*real(twopi)
          enddo
         endif
c
         do ig=1,nfilt
c
          if(dopower)then
           data(1)=pdata(1,ig,idate)       ! zero frequency
           data(2)=0.
           datar  =data(2)*sin(phi(1))+data(1)*cos(phi(1))
           datai  =data(2)*cos(phi(1))-data(1)*sin(phi(1))
           data(1)=datar
           data(2)=datai
           data(2)=0.      ! so that output is a real function
           do i=3,mxpl,2
            data(i)  =pdata(i/2+1,ig,idate)
            data(i+1)=0.
            datar    =data(i+1)*sin(phi(i/2+1))+data(i)*cos(phi(i/2+1))
            datai    =data(i+1)*cos(phi(i/2+1))-data(i)*sin(phi(i/2+1))
            data(i)  =datar
            data(i+1)=datai
            data(mxpl*2-i+2)=+data(i)      ! conjugate complex
            data(mxpl*2-i+3)=-data(i+1)
           enddo
           data(mxpl+1)=0.                 ! aliased point
           data(mxpl+2)=0.
           isign=-1
           call four1(data,mxpl,isign)
          else if(doshuffle)then
           x=ran1(idum)
           idate=int(x*ndates+1.)
           x=ran1(idum)
           ical=int(x*ncals(idate)+1.)
          endif
 
          do n=1,lastid-firstid+1
           obsvis2(n+firstid-1,ig)=modvis2(n+firstid-1,ig)
     &          +dble(gasdev(idum))*viserr(n+firstid-1,ig)
           j=int((real(hours(n+firstid-1))-ptimmin(idate))
     &                /ptimstep(idate))+1
           if(j.gt.mxpl)j=mxpl
           if(j.lt.1)j=1
           if(dopower)then
            obsvis2(n+firstid-1,ig)=obsvis2(n+firstid-1,ig)*
     &       dble(1.+data(2*j-1)/real(mxpl))
           else if(doshuffle)then
            obsvis2(n+firstid-1,ig)=visdata(n+firstid-1,ig)/
     &       dble(1.+caldata(idate,ical,ig,j))
           endif
          enddo
 
         enddo               ! Next filter
c
   40    firstid=lastid+1
         lastid=firstid-1
c
c       continue here, if date did not change
        endif
        lastid=lastid+1
        ENDDO

       endif               ! choice of method 2/1,3
c
      ENDIF                ! choice of w/out calibration errors
c.......................................
c
      if(test)then
       itest=itest+1
       mode=1
       call plotvis(modplot,modvis2,modphas,mode,device)
       if(mod(itest,3).eq.0)then
        write(outc,*)'Click on window and enter command:'
        call pgcurse(cx,cy,command)
        if(command.eq.'q')then
         call pgend
         return
        else if(command.eq.'?')then
         write(outc,*)'-q- Quit'
        endif
       endif
       goto 30
      endif
      return
      end