* DREAM = Data REduction And Merging
* Version 0 -- Ns scopes, Nt times 
* 2020 Apr Keith Horne @ St.Andrews
* 2020 May KDH @ StA - multiple wavelengths
* 2020 May KDH @ StA - data with time,lam,scope labels
* 2020 May KDH @ StA - fake data with duty cycle
* 2020 May KDH @ StA - X(t) normalise to avg 0 rms 1 mag
* 2020 May KDH @ StA - convergence plot
* 2020 May KDH @ StA - Input data from files
* 2020 May KDH @ StA - chi^2/dof on plots
* 2020 May KDH @ StA - avgrmsx
* 2020 Jun KDH @ StA - running optimal average
* 2020 Jun KDH @ StA - output files
* 2020 Jun KDH @ StA - roap
* 2020 Aug KDH @ Crail - aic, bic
* 2020 Aug KDH @ Crail - grid scan
* 2020 Aug KDH @ StA - cpu
* 2020 Aug KDH @ StA - parabolic fit to 3 best
* 2020 Aug KDH @ StA - revise convergence tests
* 2020 Aug KDH @ StA - plot lightcurve residuals in separate panel
* 2020 Sep KDH @ StA - original and revised error bars on plot
* 2020 Sep KDH @ StA - skip yellow, large errors first, then smaller
* 2020 Sep KDH @ StA - gaussian prior on scl

	module dream
* max files
	parameter( maxf = 1000 )
* max header lines
	parameter( maxh = 100 )
* max times
	parameter( maxt = 30000 )
* max wavelengths
	parameter( maxw = 100 )
* max scopes
	parameter( maxs = 50 )
* max observations
	parameter( maxo = maxt * maxw )
* max iterations
	parameter( maxi = 100000 )
* max convergence tests
	parameter( maxc = 10 )
* max grid
	parameter( maxg = 1000 )
* max plot
	parameter( maxp = maxi )

* time-series observations
	integer*4 iw_o( maxo )
	integer*4 is_o( maxo )	
	real*8   hjd_o( maxo )
	real*4     t_o( maxo )
	real*4   dat_o( maxo )
	real*4   sig_o( maxo )

	real*4   fit_o( maxo )
	real*4   dif_o( maxo )
	real*4   chi_o( maxo )

* model lightcurve
	real*4    x_o( maxo )
	real*4 sigx_o( maxo )

* calibrated data
	real*4    cal_o( maxo )
	real*4 difcal_o( maxo )
	real*4 sigcal_o( maxo )
	real*4 rmscal_o( maxo )
	real*4 sigoldcal_o( maxo )

* parameters controlled by each datum
	real*4 dofroa_o( maxo )
c	real*4 dofcal_o( maxo ) ! not yet used

* lightcurves
	real*8 hjd, hjd_t( maxt )
	real*4  to_t( maxt )
	real*4 dat_t( maxt )
	real*4 sig_t( maxt )
	real*4 fit_t( maxt )
	real*4 dif_t( maxt )

* wavelengths
	character*10 name
	character*10 name_w( maxw )
	real*4 w_w( maxw )
	real*4 dat_w( maxw )
	real*4 sig_w( maxw )
	real*4 fit_w( maxw )
	real*4 dif_w( maxw )
	integer*4 ndat_w( maxw )

* scopes
	character*10 name_s( maxs )
	integer*4 kolor_s( maxs )
	real*4 dat_s( maxs )
	real*4 sig_s( maxs )
	real*4 fit_s( maxs )
	real*4 dif_s( maxs )
	integer*4 ndat_s( maxs )

* raw data
	integer*4 nt_s( maxs )
	real*4   t_ts( maxt, maxs )
	real*4 dat_ts( maxt, maxs )
	real*4 sig_ts( maxt, maxs )

* per scope and wavelength
	real*4 dat_tsw( maxt, maxs, maxw )
	real*4 sig_tsw( maxt, maxs, maxw )
	integer*4 ndat_sw( maxs, maxw )

* calibrated data
	real*4    cal_t( maxt )
	real*4 sigcal_t( maxt )
	real*4 rmscal_t( maxt )
	real*4 sigoldcal_t( maxt )
* per scope
	real*4    cal_ts( maxt, maxs )
	real*4 difcal_ts( maxt, maxs )
	real*4 sigcal_ts( maxt, maxs )
	real*4 rmscal_ts( maxt, maxs )
	real*4 sigoldcal_ts( maxt, maxs )
* per scope and wavelength
	real*4    cal_tsw( maxt, maxs, maxw )
	real*4 difcal_tsw( maxt, maxs, maxw )
	real*4 sigcal_tsw( maxt, maxs, maxw )
	real*4 rmscal_tsw( maxt, maxs, maxw )
	real*4 sigoldcal_tsw( maxt, maxs, maxw )

* model lightcurve (running optimal average)
	real*4 t_t( maxt )
	real*4 x_t( maxt )
	real*4 sigx_t( maxt )

	integer*4 nt_w( maxw )
	real*4    t_tw( maxt, maxw )
	real*4    x_tw( maxt, maxw )
	real*4 sigx_tw( maxt, maxw )
* changes
	real*4   dx_tw( maxt, maxw )

* roa window function
	character*1 shape
	real*4 win_t( maxt )

* calibration and noise model parameters (per scope)
	real*4 add_s( maxs )
	real*4 scl_s( maxs )
	real*4 rms_s( maxs )

* per scope and wavelength
	real*4 add_sw( maxs, maxw )
	real*4 scl_sw( maxs, maxw )
	real*4 rms_sw( maxs, maxw )
* Gaussian prior
	logical priorscl
	real*4 prioravgscl_sw( maxs, maxw )
	real*4 priorrmsscl_sw( maxs, maxw )

* changes
	real*4 dadd_sw( maxs, maxw )
	real*4 dscl_sw( maxs, maxw )
	real*4 drms_sw( maxs, maxw )

* parameter error bars (per scope )
	real*4 sigadd_s( maxs )
	real*4 sigscl_s( maxs )
	real*4 sigrms_s( maxs )
* per scope and wavelength
	real*4 sigadd_sw( maxs, maxw )
	real*4 sigscl_sw( maxs, maxw )
	real*4 sigrms_sw( maxs, maxw )

* fit, residuals, normalised residuals
	real*4 fit_ts( maxt, maxs )
	real*4 dif_ts( maxt, maxs )
	real*4 chi_ts( maxt, maxs )
* per scope and wavelength
	real*4 fit_tsw( maxt, maxs, maxw )
	real*4 dif_tsw( maxt, maxs, maxw )
	real*4 chi_tsw( maxt, maxs, maxw )

* real*8 variables and arrays
	real*8 sum
	real*8 sum_t( maxt )
	real*8 sum_s( maxs )
	real*8 sum_w( maxw )

* badness of fit and degrees of freedom
	real*8 bof, chi2, bofvar, aic, bic
	real*8 dof, dofdat, doffit, dofvar, dofcal

c	real*8 bofvar_t( maxt ) ! not used
c	real*8 bofvar_s( maxs ) ! defined below
c	real*8 bofvar_w( maxw ) ! defined below

* integer arrays (e.g. for sorted values)
	integer*4 i_t( maxt )
	integer*4 i_s( maxs )
	integer*4 i_w( maxw )
	integer*4 i_o( maxo )

* -----------------------------------------------
* iteration record
	real*4 safe_i( maxi )
	real*4 tiny_i( maxi )

* per iteration, scope and wavelength
	real*4 add_isw( maxi, maxs, maxw )
	real*4 scl_isw( maxi, maxs, maxw )
	real*4 rms_isw( maxi, maxs, maxw )
	real*4 sigadd_isw( maxi, maxs, maxw )
	real*4 sigscl_isw( maxi, maxs, maxw )
	real*4 sigrms_isw( maxi, maxs, maxw )
* per iteration and case
	real*8 chi_ic( maxi, maxc )

* per iteration
	integer*4 ndat_i( maxi )
	real*4 dofcal_i( maxi )
	real*4 dofroa_i( maxi )
	real*8 bofvar_i( maxi )
	real*8   chi2_i( maxi )
	real*8    dof_i( maxi )
	real*8    bof_i( maxi )
	real*8    aic_i( maxi )
	real*8    bic_i( maxi )
	
* per wavelength
	real*4 dofcal_w( maxw )
	real*4 dofroa_w( maxw )
	real*8 bofvar_w( maxw )
	real*8   chi2_w( maxw )
	real*8    dof_w( maxw )
	real*8    aic_w( maxw )
	real*8    bic_w( maxw )
* per scope
	real*4 dofcal_s( maxs )
	real*4 dofroa_s( maxs )
	real*8 bofvar_s( maxs )
	real*8   chi2_s( maxs )
	real*8    dof_s( maxs )
	real*8    aic_s( maxs )
	real*8    bic_s( maxs )

* per iteration and wavelength
	real*4 dofcal_iw( maxi, maxw )
	real*4 dofroa_iw( maxi, maxw )
	real*8 bofvar_iw( maxi, maxw )
	real*8   chi2_iw( maxi, maxw )
	real*8    dof_iw( maxi, maxw )
	real*8    aic_iw( maxi, maxw )
	real*8    bic_iw( maxi, maxw )

* per iteration and scope
	real*4 dofcal_is( maxi, maxs )
	real*4 dofroa_is( maxi, maxs )
	real*8 bofvar_is( maxi, maxs )
	real*8   chi2_is( maxi, maxs )
	real*8    dof_is( maxi, maxs )
	real*8    aic_is( maxi, maxs )
	real*8    bic_is( maxi, maxs )

* -----------------------------------------------
* parameter grid
	real*4	grid_g( maxg )
	real*4	    hw_g( maxg )
	real*4	    dt_g( maxg )
	real*8	dofcal_g( maxg )
	real*8	dofroa_g( maxg )
	real*8	bofvar_g( maxg )
	real*8	  chi2_g( maxg )
	real*8	   bof_g( maxg )
	real*8	   dof_g( maxg )
	real*8	   aic_g( maxg )
	real*8	   bic_g( maxg )
	integer*4 nits_g( maxg )
	integer*4 ndat_g( maxg )
* calibration grid
	real*4	scl_gs( maxg, maxs )
	real*4	add_gs( maxg, maxs )
	real*4	rms_gs( maxg, maxs )
	real*4	scl_gsw( maxg, maxs, maxw )
	real*4	add_gsw( maxg, maxs, maxw )
	real*4	rms_gsw( maxg, maxs, maxw )
* calibration uncertainties
	real*4	sigscl_gs( maxg, maxs )
	real*4	sigadd_gs( maxg, maxs )
	real*4	sigrms_gs( maxg, maxs )
	real*4	sigscl_gsw( maxg, maxs, maxw )
	real*4	sigadd_gsw( maxg, maxs, maxw )
	real*4	sigrms_gsw( maxg, maxs, maxw )

* steps (convergence test names)
	character*50 name_c( maxc )
* agn
	character*20 agn
	character*20 prefix
* file names
	character*128 file, obsfile
	character*128 file_f( maxf )
	character*128 infile_sw( maxs, maxw )
* header lines
	character*128 h, h_h( maxh )
* pgplot
	real*4 plot1( maxp )
	real*4 plot2( maxp )
	real*4 plot3( maxp )
	real*4 plot4( maxp )
	real*4 plot5( maxp )
	real*4 plot6( maxp )
	real*4 plot7( maxp )
	real*4 plot8( maxp )
	real*4 plot9( maxp )
	real*4 plot10( maxp )
	real*4 plotx( maxp )
	real*4 ploty( maxp )
	real*4 plot1_s( maxp, maxs )
	real*4 plot2_s( maxp, maxs )
	real*4 plot3_s( maxp, maxs )
	character*75 xlabel, ylabel, title, label
* toggles
	logical gabby, yes, converged, ok, aok, lcnorm

	end module
*****************************************************
* global variables
	use dream
* local variables
	logical pgagain
	character*10 reply
	character*50 version

	version = '2020-Sep-09'
	write(*,*)
	write(*,*) '**************************************'
	write(*,*) '*                                    *'
	write(*,*) '* DREAM = Data REduction And Merging *'
	write(*,*) '*                                    *'
	write(*,*) '* Inter-Calibration of Lightcurves   *'
	write(*,*) '*                                    *'
	write(*,*) '* Keith Horne @ St.Andrews           *'
	write(*,*) '*   v0.0 - test with fake data       *'
	write(*,*) '*   v0.1 - multiple wavelengths      *'
	write(*,*) '*   v0.2 - general time sampling     *'
	write(*,*) '*   v0.3 - fake data duty cycle      *'
	write(*,*) '*   v0.4 - data from ascii files     *'
	write(*,*) '*   v0.5 - running optimal average   *'
	write(*,*) '*   v0.6 - excess variance           *'
	write(*,*) '*   v0.7 - output files              *'
	write(*,*) '*   v0.8 - ROA parameter count       *'
	write(*,*) '*   v0.9 - AIC and BIC               *'
	write(*,*) '*   v1.0 - Grid Scan                 *'
	write(*,*) '*   v1.1 - revised convergence test  *'
	write(*,*) '*   v1.2 - revised lightcurve plot   *'
	write(*,*) '*   v1.3 - gaussian prior on scl     *'
	write(*,*) '**************************************'
	write(*,*) version(:len1(version))

* pgplot colour indices
	kblack = 1
	kred = 2
	kgreen = 3
	kblue = 4
	kcyan = 5
	kpink = 6
	kyellow = 7
	korange = 8
	kyelgrn = 9
	kblugrn = 10
	kltblu = 11
	kpurple = 12
	kpuce = 13
	kdkgrey = 14
	kltgrey = 15
* pgplot symbols
	istar = 18
	idot = 1
	icirc = 21

	pi = 4. * atan2( 1., 1. )
	iseed = 19571095

* initialise
	no = 0
* times
	nt = 0
* scopes
	ns = 0
* wavelengths
	nw = 0
* iterations
	iter = 0
	newits = 0
	cpuiter = 0.
* priors
	priorscl = .true.
	prioravgscl = 1.d0
	priorrmsscl = 1.d-2
	priorrmsscl = 5.d-2

* lightcurve normalisation option (true if <L>=0 <L^2>=1)
	lcnorm = .true.
	lcnorm = .false.
* model timescale for running optimal average
	hw = 2.
	shape = 'C'
* model time spacing
	dt = hw / 10.

* roa halfwidth grid
	nhw = 5
	newitg = 10
	hw1 = 0.5
	hw2 = 4.0
	hwaic = 0.d0
	hwbic = 0.d0
	cpugrid = 0.

* output file prefix
	prefix = ' '
* duty (for fake data)
	duty = 0.
* gabby 
	gabby = .true.
	gabby = .false.
* convergence threshold
	tiny = 1.e-4
* safety
	safe = 0.5
	safe = 0.9

* plot symbol on lightcurves
	isymb = idot
* length of error bar terminals
	term = 0.5
	term = 1.0

	obsfile = ' '

c	goto 20

* menu
   10	write(*,*)
	write(*,*) 'DREAM menu:'
	write(*,*) ' D ... Data from file ', obsfile(:len1(obsfile))
	write(*,*) ' F ... Fake data Duty', duty
* commands requiring data
      if( no .gt. 0 ) then
	call mnmx( no, t_o, tmn, tmx )
	write(*,*) '       No', no, ' JD', tmn, tmx
	write(*,*) '       Nw', nw, ( ' ', name_w(i)(:len1(name_w(i))), i=1,nw)
	write(*,*) '       Ns', ns, ( ' ', name_s(i)(:len1(name_s(i))), i=1,ns)
* lightcurve model
	write(*,*) ' M ... Model Nt', nt, ' dt', dt
	write(*,*) ' R ... ROA shape ', shape(:len1(shape)), ' HW', hw
     &		, ' =', hw/dt, ' dt'
	write(*,*) ' N ... Normalise Lightcurve to avg 0 rms 1 ', lcnorm
* prior
      if( priorscl ) then
	write(*,*) ' P ... Prior scl', priorscl
     &		, ' avg', prioravgscl, ' rms', priorrmsscl
      else
	write(*,*) ' P ... Prior scl', priorscl
      end if

* half-width grid scan
      if( dt .gt. hw/10. ) then
	write(*,*) '** WARNING: HW / dt =', hw / dt, ' **'
      end if
	write(*,*) ' G ... Grid Scan HW', nhw, hw1, hw2

* iterations
      if( converged ) then
	write(*,*) ' I ... Iterate', iter, ' Converged :)'
      else
	write(*,*) ' I ... Iterate', iter, ' NOT Converged.'
      end if

      if( iter .gt. 0 ) then
	write(*,*) ' Z ... Zap previous iterations'
	aic = aic_i( iter )
	bic = bic_i( iter )
	chi2 = chi2_i( iter )
	bofvar = bofvar_i( iter )
	bof = chi2 + bofvar
      end if

* KDH : should add more summary diagnostics here
	write(*,*) '       Chi^2', real( chi2 ), ' bofvar', real( bofvar )
	write(*,*) '       -2lnL', real(bof), ' AIC', real(aic), ' BIC', real(bic)
	write(*,*) '  S ... Safe', safe
	write(*,*) '  T ... Tiny', tiny
* plots
	write(*,*) ' L ... Lightcurve Plot'

      if( iter .ne. 0 ) then
      if( converged ) then
	write(*,*) ' X ... Convergence Plot iterations', iter, ' Converged :)'
      else
	write(*,*) ' I ... Convergence Plot iterations', iter, ' NOT Converged.'
      end if
      end if

	write(*,*) ' C ... Calibration summary'
* output
	write(*,*) ' O ... Output to files, prefix: ' // prefix(:len1(prefix))
      end if
* fini
	write(*,*) ' Gabby diagnostics ', gabby
	write(*,*) ' Q ... Quit'

* command node
   11	call inq( 'DREAM>', ' ', reply )
	call word1( reply, l1, l2 )
	reply = reply( l1 : l2 )
	call upcase( reply, reply )

* quit
      if( reply .eq. 'Q' ) then
	write(*,*)
	write(*,*) '     Sweet DREAMs'
	write(*,*)
	stop
      end if
* gabby
      if( reply .eq. 'GABBY' ) then
	reply = 'N'
	if( gabby ) reply = 'Y'
	call inq( 'Gabby diagnostics? ', reply, reply )
	call upcase( reply, reply )
	gabby = reply .eq. 'Y'
	goto 11
      end if

	if( reply .eq. 'D' ) goto 500
	if( reply .eq. 'F' ) goto 20

* commands requiring data
      if( no .gt. 0 ) then
	if( reply .eq. 'M' ) goto 30
	if( reply .eq. 'R' ) goto 40
	if( reply .eq. 'I' ) goto 50
	if( reply .eq. 'L' ) goto 100
	if( reply .eq. 'X' ) goto 200
	if( reply .eq. 'C' ) goto 300
	if( reply .eq. 'O' ) goto 400
	if( reply .eq. 'G' ) goto 600

* toggle priors
      if( reply .eq. 'P' ) then
	reply = 'N'
	if( priorscl ) reply = 'Y'
	call inq( 'Gaussian Prior on scl ', reply, reply )
	call upcase( reply, reply )
	priorscl = reply .eq. 'Y'
      if( priorscl ) then
	avg = prioravgscl
	rms = priorrmsscl
	call inq2r( ' Prior RMS and AVG on SCL :', rms, avg )
      if( avg .le. 0. .or. rms .le. 0. ) then
	write(*,*) '** ERROR : INVALID RMS or AVG', rms, avg
      else
      if( priorrmsscl .ne. rms ) then
	write(*,*) ' RMS(SCL)', priorrmsscl, ' ->', rms
	priorrmsscl = rms
      end if
      if( prioravgscl .ne. avg ) then
	prioravgscl = avg
	write(*,*) ' AVG(SCL)', prioravgscl, ' ->', avg
      end if
      end if
      end if
	goto 11
      end if

* reset iteration counter
      if( reply .eq. 'Z' ) then
	i = iter
	write(*,*) 'Iterations', iter, ' of max', maxi
	call inqi( 'Reset iterations to', iter, i )
	if( i .ge. 0 ) iter = max( 0, min( maxi, i ) )
	goto 11
      end if


* normalisation
      if( reply .eq. 'N' ) then
	reply = 'N'
	if( lcnorm ) reply = 'Y'
	call inq( 'Normalise Lightcurve to avg 0 rms 1 ?', reply, reply )
	call upcase( reply, reply )
	lcnorm = reply(1:1) .eq. 'Y'
	goto 11
      end if

* safe
      if( reply .eq. 'S' ) then
	write(*,*) 'Shorten steps by safety factor', safe
	call inqr( 'Safe', safe, val )
	if( val .gt. 0. ) safe = min( 1., val )
	goto 11
      end if
* tiny
      if( reply .eq. 'T' ) then
	write(*,*) 'Converge when step/sig < ', tiny
	call inqr( 'Tiny', tiny, val )
	if( val .gt. 0. ) tiny = min( 1., val )
	goto 11
      end if

* end if no>0
      end if


	write(*,*) 'Command not recognised. Please try again.'
	goto 10

* -------------------------------------------
* fake data for testing
   20	continue

* KDH - fake data fully fills nt times, nw wavelengths, ns scopes
* wavelengths
   	nw = 3
	w_w( 1 ) = 3000.
	w_w( 2 ) = 5000.
	w_w( 3 ) = 7000.
	name_w( 1 ) = 'u'
	name_w( 2 ) = 'g'
	name_w( 3 ) = 'r'
      do iw = 1, nw
* scopes
	ns = 5
	snr = 10.
      do is = 1, ns
	name_s( is ) = ' '
	kolor_s( is ) = is
	call append_i( name_s( is ), is )
	add_sw( is, iw ) = is + iw
	scl_sw( is, iw ) = is + iw
	rms_sw( is, iw ) = is / snr
      end do
* times
	nt = 100
	p = nt / ( 5 + iw )
	b = iw 
	a = 0.2
      do it = 1, nt
	t = it
	x = b + a * sin( 2. * pi * t / p )
	t_t( it ) = t
	t_tw( it, iw ) = t
	x_tw( it, iw ) = x
      do is = 1, ns 
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )

* fake data have rms = sig
	sig = rms
* old error bar
	sigold = sig
* new error bar
	sig = sqrt( sig * sig + rms * rms )

	fit = add + scl * x
	dat = rang( fit, sig, iseed )
	dif = dat - fit

	cal = ( dat - add ) / scl
	sigcal = sig / scl
	rmscal = rms / scl
	sigoldcal = sigold / scl
* stow
	dat_ts( it, is ) = dat
	sig_ts( it, is ) = sig
	dif_ts( it, is ) = dif

	cal_ts( it, is ) = cal
	sigcal_ts( it, is ) = sigcal
	rmscal_ts( it, is ) = rmscal
	sigoldcal_ts( it, is ) = sigoldcal

	dat_tsw( it, is, iw ) = dat_ts( it, is )
	sig_tsw( it, is, iw ) = sig_ts( it, is )   
	dif_tsw( it, is, iw ) = dif_ts( it, is )

	cal_tsw( it, is, iw ) = cal_ts( it, is )
	sigcal_tsw( it, is, iw ) = sigcal_ts( it, is )
	rmscal_tsw( it, is, iw ) = rmscal_ts( it, is )
	sigoldcal_tsw( it, is, iw ) = sigoldcal_ts( it, is )
      end do
* next time
      end do
	nt_w( iw ) = nt
* next wavelength
      end do

* report
	write(*,*) 'Fake data ----------------------------'
	write(*,*) 'Scopes:', ns, ' Times:', nt

      do iw = 1, nw
	nt = nt_w( iw )
	lam = nint( w_w( iw ) )
	name = name_w( iw )
      do is = 1, ns
	call avgrms( nt, dat_tsw( 1, is, iw ), avg, rms )
	write(*,*) 'Scope', is, ' ', name(:len1(name)), lam, avg, rms
      end do
	write(*,*) '-----------------------------'
      end do

* KDH - dilute the data
* duty cycle
	duty = 1. / ns
	write(*,*) 'Scopes', ns, ' duty', duty
	call inqr( 'Duty Cycle', duty, duty )
	if( duty .le. 0. ) goto 11

	write(*,*) '-----------------------------'

* stow fake data in the observation data bank
	n = 0
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	name = name_w( iw )
	kw = 0
      do is = 1, ns
	ks = 0
      do it = 1, nt
* gather
	t = t_t( it )
	dat = dat_tsw( it, is, iw )
	sig = sig_tsw( it, is, iw )
	dif = dif_tsw( it, is, iw )

	cal = cal_tsw( it, is, iw )
	sigcal = sigcal_tsw( it, is, iw )
	rmscal = rmscal_tsw( it, is, iw )
	sigoldcal = sigoldcal_tsw( it, is, iw )
* random duty cycle
      if( ranu( 0., 1., iseed ) .le. duty ) then
* stow
	ks = ks + 1
	 n =  n + 1
	 iw_o( n ) = iw
	 is_o( n ) = is
	  t_o( n ) = t
	dat_o( n ) = dat
	sig_o( n ) = sig
	dif_o( n ) = dif

	cal_o( n ) = cal
	sigcal_o( n ) = sigcal
	rmscal_o( n ) = rmscal
	sigoldcal_o( n ) = sigoldcal
      end if
* next time it
      end do
	kw = kw + ks
* report
	p = real( ks ) / nt
	write(*,*) 'Scope', is, ' Epochs', ks, ' duty', p
	n1 = n - ks + 1
	call avgrms( ks, dat_o( n1 ), avg, rms)
	write(*,*) name(:len1(name)), lam, ' avg,rms', avg, rms
* next scope is
      end do
* report
	p = real( kw ) / nt
	write(*,*) name(:len1(name)), lam, ' Epochs', kw, ' duty', p
	write(*,*) '-----------------------------'
* next wavelength iw
      end do
	no = n

	iter = 0
	obsfile = 'FAKE'
	converged = .false.

* KDH - do later?
	dofroa = 0.
	dofcal = 0.

* do BoF analysis
	goto 70
c	goto 11

* -------------------------------------------
   30	write(*,*) 'Reset calibration parameters'

	iter = 0

* data time span and interval
	call mnmx( no, t_o, t1, t2 )
	span = t2 - t1
	dtavg = span / max( 1, no - 1 )
* report
	write(*,*) 'Observation times:', no, t1, t2
	write(*,*) ' Mean interval:', dtavg

* expand model time span
	add = span * 0.05
	tmn = t1 - add
	tmx = t2 + add
* model time interval
	call inqr( 'Model time interval', dt, val )
	if( val .eq. 0. ) val = dtavg
	if( val .le. 0. ) goto 11
	dt = val

* even time grid
	tspan = tmx - tmn
      if( tspan .gt. 10. ) then
	tmn = int( tmn )
	tmx = int( tmx ) + 1.
	tspan = tmx - tmn
      end if
	nt = 1 + tspan / dt
	if( nt .gt. maxt ) write(*,*) '** Model times', nt, ' max', maxt
	nt = min( nt, maxt )
	dt = tspan / max( 1, nt - 1 )
	tmx = tmn + dt * max( 1, nt - 1 )
	write(*,*) 'Model times:', nt, tmn, tmx, dt

* stow model times (per filter)
      do it = 1, nt
	p = ( it - 1. ) / max( 1, nt - 1 )
	t = tmn * ( 1. - p ) + p * tmx
	t_t( it ) = t
      end do
      do iw = 1, nw
	nt_w( iw ) = nt
      do it = 1, nt
	t_tw( it, iw ) = t_t( it )
      end do
      end do

* initial add and scale parameters
	add = 0.
	scl = 1.
* wavelengths
      do iw = 1, nw
* gather data for wavelength iw
	n = 0
      do io = 1, no
      if( iw .eq. iw_o( io ) ) then
	n = n + 1
	plot1( n ) = dat_o( io )
	plot2( n ) = sig_o( io )
      end if
      end do

* avg and rms
	call avgrms( n, plot1, avg, rms )
	call avgrmsx( n, plot1, plot2, good, avgx, rmsx, sigavgx, sigrmsx )

	name = name_w( iw )
	write(*,*) '--------------- ', name(:len1(name)), ' --------------'
	write(*,*) ' data', n, ' avg', avg, ' rms', rms
	write(*,*) ' good', good, ' avgx', avgx, ' rmsx', rmsx

	avg = avgx
	rms = rmsx

* reset model lightcurve
	if( n .lt. 1 ) avg = 1. 
c	if( n .lt. 2 ) rms = 1.
c	rms = 0.
      do it = 1, nt
	t_tw( it, iw ) = t_t( it )
	x_tw( it, iw ) = avg
      end do
	write(*,*) 'Lightcurve model set to', avg
	dofroa = 1.

* calibration parameters
      do is = 1, ns
	add_sw( is, iw ) = add
	scl_sw( is, iw ) = scl
	rms_sw( is, iw ) = rms
	name = name_s( is )
	write(*,*) ' Scope ', name(:len1(name)), ' add', add
	write(*,*) ' Scope ', name(:len1(name)), ' scl', scl
	write(*,*) ' Scope ', name(:len1(name)), ' rms', rms
      end do
* next wavelength
      end do
	write(*,*) 'Calibration Parameters Reset.'
	dofcal = 0.
* stow
	nt_w( iw ) = nt
	dofroa_w( iw ) = dofroa
	dofcal_w( iw ) = dofcal

* report
	write(*,*) 'Model times:', nt, tmn, tmx
	write(*,*) 'Model interval:', dt

* ROA shape
   40	write(*,*) 'Running Optimal Average Window Shape'
	write(*,*) ' L ... Lorentzian W(x) = 1/(1+x^2)'
	write(*,*) ' G ... Gaussuan   W(x) = exp(-x^2/2)'
	write(*,*) ' B ... Boxcar     W(x) = 1 for |x|<1'
	write(*,*) ' C ... 1/Cosh     W(x) = 1/cosh(x)'

	call inq( 'ROA Window Shape ', shape, reply )

	call word1( reply, l1, l2 )
	reply = reply( l1 : l2 )
	call upcase( reply, reply )

      if( reply .eq. 'L' ) then
	shape = 'L'
      else if( reply .eq. 'G' ) then
	shape = 'G'
      else if( reply .eq. 'B' ) then
	shape = 'B'
      else if( reply .eq. 'C' ) then
	shape = 'C'
      else
	write(*,*) '** Sorry, not Recognised. Please try Again.'
	goto 11
      end if

* ROA halfwidth
	val = hw
	call inqr( 'ROA Half-Width', hw, val )
	if( val .le. 0. ) goto 11
	hw = val

	if( newits .le. 0 ) goto 11
	goto 50

* ------------------------------------------------
* grid scan
  600	write(*,*) '----------------------------------------'
	write(*,*) 'DREAM Grid Scan Menu:'
	write(*,*) ' H ... Half-Width grid:', nhw, hw1, hw2
* iterations
c	n = newitg
c	if( n .le. 0 ) n = newits
      if( converged ) then
	write(*,*) ' I ... Iterations done', itg, ' Converged :)'
      else
	write(*,*) ' I ... Iterations done', itg, ' NOT converged.'
      end if

      if( hwaic .gt. 0.d0 ) then
	write(*,*) ' A ... Center on best AIC', real(bestaic), ' HW', real(hwaic)
      end if
      if( hwbic .gt. 0.d0 ) then
	write(*,*) ' B ... Center on best BIC', real(bestbic), ' HW', real(hwbic)
      end if
	write(*,*) ' C ... Center on current HW', real(hw)
      if( nhw .gt. 1 ) then
	if( xpand .le. 0. ) xpand = 0.1
	write(*,*) ' X ... Expand grid by factor', xpand
      end if
c	write(*,*) ' G ... Grid Scan (', nhw, ' x', n, ') =', nhw * n
      if( itg .gt. 0 ) then
	write(*,*) ' P ... Plot Grid Scan'
      end if
	write(*,*) ' Q ... Quit Grid Scan'

  610	call inq( 'DREAM Grid Scan > ', ' ', reply )

	call upcase( reply, reply )
	call word1( reply, l1, l2 )
	reply = reply( l1 : l2 )

	if( reply .eq. 'Q' ) goto 10
* plot
	if( reply .eq. 'P' ) goto 650

* iterate
      if( reply .eq. 'I' ) then
	write(*,*) 'Iterations done', itg

	call inqi( 'New iterations', newitg, i )
	if( i .le. 0 ) goto 610
	newitg = i
	goto 51

* grid range
      else if( reply .eq. 'H' ) then
	xmn = hw1
	xmx = hw2
	call inq2r( 'Grid Halfwidth range', xmn, xmx )
      if( xmn .gt. 0.d0 .and. xmx .gt. xmn ) then
	hw1 = min( xmn, xmx )
	hw2 = max( xmn, xmx )
      end if

* number of grid points
	i = nhw
	write(*,*) 'Max Grid points', maxg
	write(*,*) 'Recommend 3 5 9 17 33'
	call inqi( 'Number of Grid points', i, i )
	if( i .gt. 0 ) nhw = min( maxg, i )
	itg = 0

* expand grid
      else if( reply .eq. 'X' ) then
	if( xpand .le. 0.d0 ) xpand = 1.d0
	x = xpand
	call inqr( 'Expand grid by factor', x, x )
	if( x .le. 0. ) goto 600
	div = exp( abs( alog( hw2 / hw1 ) ) * x / 2. )
	cen = sqrt( hw1 * hw2 )
	hw1 = cen / div
	hw2 = cen * div
	itg = 0

* recentre the grid
      else if( reply .eq. 'A'
     &	  .or. reply .eq. 'B'
     &	  .or. reply .eq. 'C' ) then
	cen = hw
	if( reply .eq. 'A' ) cen = hwaic
	if( reply .eq. 'B' ) cen = hwbic
	if( cen .le. 0. ) cdn = hw
	div = exp( abs( alog( hw2 / hw1 ) ) / 2. )
	hw1 = cen / div
	hw2 = cen * div
	itg = 0

* Grid scan ----------------
c      else if( reply .eq. 'G' .or. reply .eq. 'M' ) then
* new parameter grid 
c       if( itg .le. 0 .or. reply .ne. 'M' ) then
c	goto 51
c      end if
* now run it
c	goto 52

* command not recognised
      else
	write(*,*) 'Sorry, command not recognised. Please try again.'
      end if
	goto 600

* -------------------------------------------
   50	write(*,*) 'Iterations done', iter
	call inqi( 'New iterations:', newits, newits )
	if( newits .le. 0 ) goto 11

* define 1-point grid scan
	ngrid = 1
	grid_g( 1 ) = hw
	last = newits
	goto 52

* set up parameter grid 
   51	g1 = alog( hw1 )
	g2 = alog( hw2 )
	ngrid = nhw
	last = min( maxi, newitg )
      do i = 1, ngrid
	p = ( i - 1.d0 ) / max( 1, nhw-1 )
	g = exp( g1 * ( 1. - p ) + p * g2 )
	hw_g( i ) = g
	grid_g( i ) = g
      end do
* reset iteration
	last = newitg
	goto 52

* -------------------------------------------
* run the grid scan
   52	write(*,*) 'Grid Scan HW', ngrid, grid_g(1), grid_g(ngrid)
	write(*,*) 'Previous iterations', itg
	write(*,*) 'New iterations per grid point', last
	itg0 = itg
	cpudone = 0.
* grid scan
      do ig = 1, ngrid
	grid = grid_g( ig )
* halfwidth
	hw = grid
* grid iteration counter
	itg = itg0
* cpu counter
	cputime = 0.
	call cpu( dcpu )
* iterations
      do its = 1, last

* increment iteration counter
      if( ngrid .gt. 1 ) then
	itg = min( maxi, itg + 1 )
	iti = itg
      else
	iter = min( maxi, iter + 1 )
	iti = iter
      end if

* define algorithm
* dimensionless lightcurve, normalised to <X>=0 and <X^2>=1
      if( lcnorm ) then
	isteplc  = 1
	normlc   = 1
	istepadd = 2
	istepscl = 3
	isteprms = 4
	normadd = 0
	normscl = 0
	nstep = 4

* dimensionless scale factors, <scl>=1 and <add>=0
      else
* lightcurve
	isteplc  = 1
	normlc  = 0
* add ( norm to 0 )
	istepadd = 2
	normadd  = istepadd
* scl ( normal to 1 )
	istepscl = 3
	normscl = istepscl
c	if( priorscl ) normscl = 0
* rms
	isteprms = 4
	nstep = 4
      end if

* iteration reports
	yes = gabby .or. its .eq. 1
* converged until otherwise
	converged = .true.

	if( yes ) write(*,*) 'New iteration', its, ' of', last

* do the steps -----------------------------------------
      do istep = 1, nstep

* evaluate current model and residuals
      do io = 1, no
	  t =   t_o( io )
	dat = dat_o( io )
	sig = sig_o( io )
	 iw =  iw_o( io )
	 is =  is_o( io )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
* old error bar
	sigold = sig
* new error bar 
	sig = sqrt( sig * sig + rms * rms )
* interpolate model lightcurve
	nt = nt_w( iw )
	call lookup( nt, t_tw( 1, iw ), t, i1, i2, p )
	x = terp1(   nt, x_tw( 1, iw ),    i1, i2, p )
	fit = add + scl * x
	cal = ( dat - add ) / scl
	dif = dat - fit

* KDH : RMS IS FOR RAW DATA, RMSCAL FOR CALIBRATED
	difcal = cal - x
	sigcal = sig / scl
	rmscal = rms / scl
	sigoldcal = sigold / scl
* stow
	  x_o( io ) = x
	fit_o( io ) = fit
	dif_o( io ) = dif

	cal_o( io ) = cal
	difcal_o( io ) = difcal 
	sigcal_o( io ) = sigcal
	rmscal_o( io ) = rmscal
	sigoldcal_o( io ) = sigoldcal
* next obsn io
      end do

* -------------------------------------------
* update noise model : rms_sw( is, iw ) for ns scopes nw filters
      if( istep .eq. isteprms ) then	
	if( yes ) write(*,*) 'Revise RMS'
	name_c( istep ) = 'rms'
	aok = .true.
	chi = 0.
      do iw = 1, nw
	if( yes ) write(*,*) 'Wavelength', iw, w_w( iw )
      do is = 1, ns
* gather
	n = 0
      do io = 1, no
      if( iw_o( io ) .eq. iw .and. is_o( io ) .eq. is ) then
	n = n + 1
	dif_t( n ) = dif_o( io )
	sig_t( n ) = sig_o( io )
      end if
      end do
	if( yes ) write(*,*) 'Scope', is, ' data', n
* no data
      if( n .le. 0 ) then
	write(*,*) '** No relevant data for scope', is
      else
* revise
	old = rms_sw( is, iw )
* rms of extra variations
	call avgrmsx( n, dif_t, sig_t, good, avg, rms, sigavg, sigrms )
      if( gabby ) then
	write(*,*) 'avgrmsx(', n, good, avg, rms, sigavg, sigrms, ' )'
      end if
	change = rms - old
	change = change * safe
	rms = old + change
c	sig = rms * sqrt( 2. / max( 1, n - 1 ) )
	sig = sigrms
* stow
	   rms_sw( is, iw ) = rms
	sigrms_sw( is, iw ) = sig
* converged?
	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
	lam = nint( w_w( iw ) )
      if( gabby ) then
	write(*,*) 'rms(', is, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
      end if
* next scope is
      end do
* next wavelength iw
      end do
	converged = converged .and. aok
* stow
	chi_ic( iti, istep ) = chi

      end if
* -------------------------------------------
* revise add_s( is ) for ns scopes
      if( istep .eq. istepadd ) then
	if( yes ) write(*,*) 'Revise ADD'
	name_c( istep ) = 'add'
	aok = .true.
	chi = 0.
* bands
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	if( yes ) write(*,*) 'Wavelength', iw, w_w( iw )
* scopes
      do is = 1, ns
	rms = rms_sw( is, iw )
* gather relevant data
	n = 0
      do io = 1, no
      if( iw_o( io ) .eq. iw .and. is_o( io ) .eq. is ) then
	dif = dif_o( io )
	sig = sig_o( io )
	sig = sqrt( sig * sig + rms * rms )
	n = n + 1
	dif_t( n ) = dif
	sig_t( n ) = sig
      end if
      end do
	if( yes ) write(*,*) 'Scope', is, ' data', n
      if( n .le. 0 ) then
      else
* revise
	old = add_sw( is, iw )
	call optavg( n, dif_t, sig_t, change, sig )
	change = change * safe
	add = old + change
* stow
	   add_sw( is, iw ) = add
	sigadd_sw( is, iw ) = sig
	  dadd_sw( is, iw ) = change

* test convergence

	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'add(', is, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
      end if

* next scope is
      end do
* next band iw
      end do

* convergence test if no normadd
      if( istep .ne. normadd ) then
	converged = converged .and. aok
      end if
	chi_ic( iti, istep ) = chi

      end if
* -------------------------------------------
*  normalise add_s to mean of 0

       if( istep .eq. normadd ) then

	if( yes ) write(*,*) 'Normalise Offsets <ADD>=0'
	name_c( istep ) = '<add>'

	chi = 0.
	aok = .true.
* bands
      do iw = 1, nw
	lam = nint( w_w( iw ) )
c	call avgrms( ns, add_sw(1,iw), avg, rms )
	call optavg( ns, add_sw(1,iw), sigadd_sw(1,iw), avg, sig )
* scopes
      do is = 1, ns
* reap
	old =    add_sw( is, iw )
	sig = sigadd_sw( is, iw )
* update
	change = - avg
	change = change * safe
	add = old + change

* combine changes from add and <add>
	oldchange = dadd_sw( is, iw )
	change = change + oldchange
	
* stow
	add_sw( is, iw ) = add
	dadd_sw( is, iw ) = change

* converged?
	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'add(', is, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
* next scope
      end do
* next band
      end do
* 
	converged = converged .and. aok
	chi_ic( iti, istep ) = chi

      end if
* -------------------------------------------
* revise scl_s( is ) for ns scopes
      if( istep .eq. istepscl ) then
	if( yes ) write(*,*) 'Revise SCL'
	name_c( istep ) = 'scl'
	chi = 0.d0
	aok = .true.
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	if( yes ) write(*,*) 'Wavelength', iw, w_w( iw )
      do is = 1, ns
* gather relevant data
	n = 0
      do io = 1, no
      if( iw_o( io ) .eq. iw .and. is_o( io ) .eq. is ) then
	t = t_o( io )
	nt = nt_w( iw )
	call lookup( nt, t_tw(1,iw), t, i1, i2, p )
	n = n + 1
	  x_t( n ) = terp1( nt, x_tw(1,iw), i1, i2, p )
	dif_t( n ) = dif_o( io )

	rms = rms_sw( is, iw )
	sig = sig_o( io )
	sig = sqrt( sig * sig + rms * rms )

	sig_t( n ) = sig
      end if
      end do
	if( yes ) write(*,*) 'Scope', is, ' wave', iw, ' data', n

* gaussian prior on scl 
      if( priorscl ) then
	prioravg = prioravgscl_sw( is, iw )
	priorrms = priorrmsscl_sw( is, iw )
* KDH hardwire for testing
	prioravg = prioravgscl
	priorrms = priorrmsscl

* old pre-prior result
      if( yes ) then
	datold =    scl_sw( is, iw )
	sigold = sigscl_sw( is, iw )
	chiold = ( datold - prioravg ) / priorrms
* new pre-prior result
	call optscl( n, dif_t, sig_t, x_t, change, signew )
	datnew = datold + change
	chinew = ( datnew - prioravg ) / priorrms
* report
	write(*,*) 'Scope', is, ' wave', iw, ' data', n
	write(*,*) 'Prior scl', prioravg, ' +/-', priorrms
	write(*,*) '  Old scl', datold, ' +/-', sigold, ' chi', chiold
	write(*,*) '  New scl', datnew, ' +/-', signew, ' chi', chinew
      end if

* gaussian prior acts as 1 more datum for optimal scaling
	scl = scl_sw( is, iw )
	n = n + 1
	x_t( n ) = 1.d0
	dif_t( n ) = prioravg - scl
	sig_t( n ) = priorrms
      end if

      if( n .le. 0 ) then
      else
* revise
	old = scl_sw( is, iw )
	call optscl( n, dif_t, sig_t, x_t, change, sig )
	change = change * safe
	scl = old + change

      if( priorscl ) then
      if( yes ) then
	chiwithprior = ( scl - prioravg ) / priorrms
	write(*,*) '  --> scl', scl, ' +/-', sig, ' chi', chiwithprior
      end if
      end if

* keep positive
      if( scl .le. 0. ) then 
	write(*,*) '** WARNING : STEP REDUCED TO KEEP SCL >0'
	scl = abs( old ) / 2. 
	change = scl - old
      end if
* stow
	  dscl_sw( is, iw ) = change
	   scl_sw( is, iw ) = scl
	sigscl_sw( is, iw ) = sig

* test convergence
	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'scl(', is, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
      end if
* next scope is
      end do
* next wavelength iw
      end do

* convergence test if no normscl
      if( istep .ne. normscl ) then
	converged = converged .and. aok
      end if
	chi_ic( iti, istep ) = chi
      end if
* -------------------------------------------
* normalise scl_s to mean of 1
      if( istep .eq. normscl ) then
	if( yes ) write(*,*) 'Normalise Scale Factors <SCL> = 1'
	name_c( istep ) = '<scl>'
	aok = .true.
	chi = 0.
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	call avgrms( ns, scl_sw(1,iw), avg, rms )
      do is = 1, ns
* reap
	old =    scl_sw( is, iw )
	sig = sigscl_sw( is, iw )
* update
	scl = old / avg
	change = scl - old
	change = change * safe
	scl = old + change
* stow
	scl_sw( is, iw ) = scl

* combine changes from scl and <scl>
	oldchange = dscl_sw( is, iw )
	change = change + oldchange
	dscl_sw( is, iw ) = change

* converged?
	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'scl(', is, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
      end do
      end do
	converged = converged .and. aok
	chi_ic( iti, istep ) = chi
      end if
* -------------------------------------------
* revise model lightcurve
      if( istep .eq. isteplc ) then
	if( yes ) write(*,*) 'Revise Lightcurve X(t)'
	name_c( istep ) = 'lcv'
	chi = 0.
	aok = .true.
* wavelengths
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	name = name_w( iw )
	nt = nt_w( iw )
	if( yes ) write(*,*) 'Wavelength', iw, lam, name(:len1(name)), ' times', nt

* gather data from all scopes for wavelength iw
	n = 0 
	lost = 0
      do io = 1, no
      if( iw_o( io ) .eq. iw ) then
      if( n .lt. maxt ) then
* gather
	is = is_o( io )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
	  t =   t_o( io )
	dat = dat_o( io )
	sig = sig_o( io )
* calculate
	sig = sqrt( sig * sig + rms * rms )
	x = ( dat - add ) / scl
	sigx = sig / scl
* stow
	n = n + 1
	  to_t( n ) = t
	   x_t( n ) = x
	sigx_t( n ) = sigx
	   i_t( n ) = io
      else
	lost = lost + 1
      end if
      end if
* next io
      end do	
	if( yes ) write(*,*) 'Data', n, ' lost', lost
	if( lost .gt. 0 ) stop
* time loop it
      do it = 1, nt
	old = x_tw( it, iw )
	t = t_tw( it, iw )
* running optimal average at time t
	call roa( shape, hw, n, to_t, x_t, sigx_t, 1, t, avg, sig )
* revise x
	change = avg - old
	change = change * safe
	x = old + change
* stow
	   x_tw( it, iw ) = x
	sigx_tw( it, iw ) = sig
	  dx_tw( it, iw ) = change
* test convergence
	ok = ok .and. abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'x(', it, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
* next time it
      end do

* roa lightcurve parameter count
	if( iw .eq. 1 ) dofroa = 0.
	sum = 0.d0
      do i = 1, n
	io = i_t( i )
	dofroa = roap1( shape, hw, n, to_t, sigx_t, io )
	dofroa_o( io ) = dofroa
	sum = sum + dofroa
      end do
	dofroa_w( iw ) = sum
	dofroa = dofroa + sum
      if( gabby ) then
	write(*,*) 'ROA lightcurve parameters', real( sum )
      end if
	
* next wavelength iw
      end do

* convergence test if no normlc
      if( istep .ne. normlc ) then
	converged = converged .and. aok
      end if
	chi_ic( iti, istep ) = chi
      end if
* -------------------------------------------
* normalise lightcurve to mean 0 rms 1
      if( istep .eq. normlc ) then
	if( yes ) write(*,*) 'Normalise Lightcurve to mean 0 rms 1 mag'
	name_c( istep ) = '<lcv>'
	aok = .true.
	chi = 0.
* wavelenths
      do iw = 1, nw
	lam = nint( w_w( iw ) )
	name = name_w( iw )
	if( yes ) write(*,*) '-----------------------', lam, ' ', name(:len1(name))
* skip no data
	nt = nt_w( iw )
      if( nt .lt. 0 ) then
	write(*,*) '** No data', nt
	exit
      end if
* report
	call mnmx( nt, x_tw( 1, iw ), xmn, xmx )
	call avgrms( nt, x_tw( 1, iw ), avg, rms )
	write(*,*) ' X(t) min,max', xmn, xmx
	write(*,*) '      avg,rms', avg, rms
* gather
      do it = 1, nt
	plot1( it ) =    x_tw( it, iw )
	plot2( it ) = sigx_tw( it, iw )
	plot3( it ) =   dx_tw( it, iw )
      end do
* avg and rms
	call mnmx( nt, plot1, xmn, xmx )
	call avgrms( nt, plot1, avg, rms )
	call avgrms( nt, plot3, davg, drms )

      if( yes ) then
	write(*,*) '  X(t) min,max', xmn, xmx
	write(*,*) '       avg,rms', avg, rms
	write(*,*) '     davg,drms', davg, drms
      end if

* normalise to avg = 0 and rms = 1
	call avgrms( nt, plot1, avg, rms )
* times
      do it = 1, nt
* reap
	old = plot1( it )
	sig = plot2( it )
* update
	x = old - avg
      if( rms .gt. 0. ) then
	x = x / rms
	sig = sig / rms
      end if
	change = x - old
	change = change * safe
	x = old + change	

	oldchange = dx_tw( it, iw )
	change = change + oldchange
	dx_tw( it, iw ) = change

* stow
	plot1( it ) = x
	plot2( it ) = sig


* converged?
	ok = abs( change ) .lt. sig * tiny
	aok = aok .and. ok
      if( gabby ) then
	write(*,*) 'X(', it, lam, ' )=(', old, ' +', change, ' ) +/-', sig, ok
      end if
	chi = max( chi, abs( change / sig ) )
* next time it
      end do
* stow
      do it = 1, nt
	   x_tw( it, iw ) = plot1( it )
	sigx_tw( it, iw ) = plot2( it )
      end do
* next wavelength iw
      end do	

	converged = converged .and. aok
	chi_ic( iti, istep ) = chi
      end if

* report
	name = name_c( istep )
	write(*,*) 'Iter', iti, ' Step ', chi, ' ', name(:len1(name))

* ------------------------------------------
* calibrated data, current model and residuals
      do io = 1, no
	  t =   t_o( io )
	dat = dat_o( io )
	sig = sig_o( io )
	 iw =  iw_o( io )
	 is =  is_o( io )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
* old error bar
	sigold = sig
* new error bar
	sig = sqrt( sig * sig + rms * rms )
* interpolate model lightcurve
	nt = nt_w( iw )
	call lookup( nt, t_tw( 1, iw ), t, i1, i2, p )
	x = terp1(   nt, x_tw( 1, iw ),    i1, i2, p )
	fit = add + scl * x
	dif = dat - fit

	cal = ( dat - add ) / scl
	difcal = cal - x
	sigcal = sig / scl
	rmscal = rms / scl
	sigoldcal = sigold / scl
* stow
	  x_o( io ) = x
	fit_o( io ) = fit
	dif_o( io ) = dif

	cal_o( io ) = cal
	difcal_o( io ) = difcal 
	sigcal_o( io ) = sigcal
	rmscal_o( io ) = rmscal
	sigoldcal_o( io ) = sigoldcal 
* next obsn io
      end do

* next step
      end do

* BoF analysis -------------------------------------

* zero sums and counters
   70	ndat = 0
	sum = 0.d0
	bofvar = 0.d0
	sumdofroa = 0.d0
* per wavelength
      do iw = 1, nw
	sum_w( iw ) = 0.d0
	ndat_w( iw ) = 0
	dofroa_w( iw ) = 0.d0
	bofvar_w( iw ) = 0.d0
      end do
* per scope
      do is = 1, ns
	sum_s( is ) = 0.d0
	ndat_s( is ) = 0
	dofroa_s( iw ) = 0.d00
	bofvar_s( is ) = 0.d0
* per scope and wavelength
      do iw = 1, nw
	ndat_sw( is, iw ) = 0
      end do
      end do

* data
      do io = 1, no
* skip if error bar < 0
	sig = sig_o( io )
      if( sig .gt. 0. ) then
* recall
	iw  =  iw_o( io )
	is  =  is_o( io )
	dif = dif_o( io )
	rms = rms_sw( is, iw )
	dofroa = dofroa_o( io )
* calculate
	ss = sig * sig
	rr = rms * rms
	v = ss + rr
	sig = sqrt( v )
	chi = dif / sig
* chi2 sum
	add = chi * chi
	sum = sum + add
	sum_w( iw ) = sum_w( iw ) + add
	sum_s( is ) = sum_s( is ) + add
* ln(var) sum
	add = dlog( 1.d0 + dble( rr / ss ) )
	bofvar         = bofvar         + add
	bofvar_w( iw ) = bofvar_w( iw ) + add
	bofvar_s( is ) = bofvar_s( is ) + add
* roa parameter count
	add = dofroa
	sumdofroa = sumdofroa + add
	dofroa_w( iw ) = dofroa_w( iw ) + add
	dofroa_s( is ) = dofroa_s( is ) + add
* ndat sum
	ndat = ndat + 1
	ndat_w( iw ) = ndat_w( iw ) + 1
	ndat_s( is ) = ndat_s( is ) + 1
	ndat_sw( is, iw ) = ndat_sw( is, iw ) + 1
      end if
* next datum io
      end do

	write(*,*) 'Grid', ig, ' HW', hw, ' Iter', iti, ' ---------', converged
	write(*,*) 'Data', ndat, ' of', no, ' lost', no - ndat

* calibration parameters (3 per wavelength per scope)
* assumes all scopes observed at all wavelengths

	dofcal = 3 * ns * nw
	if( iti .le. 0 ) dofcal = 0
* gather
	chi2   = sum
	dofroa = sumdofroa
* compute
	doffit = dofroa + dofcal
	dof = ndat - doffit
	bof = chi2 + bofvar
	aic = bof + doffit * 2.
	bic = bof + doffit * alog( real( ndat ) )	
* report
	write(*,*) 'Chi2', chi2, ' /', dof, ' =', chi2/dof
	write(*,*) ' Sum(ln(var))', bofvar, ' RoA params', dofroa
	write(*,*) '  AIC', aic, '  BIC', bic
	i = iti - 1
      if( i .gt. 0 ) then
	write(*,*) ' dAIC', aic-aic_i(i), ' dBIC', bic-bic_i(i)
      end if

* stow per-iteration results
	i = iti
      if( i .gt. 0 ) then
	  safe_i( i ) = safe
	  tiny_i( i ) = tiny
	dofcal_i( i ) = dofcal
	dofroa_i( i ) = dofroa
	bofvar_i( i ) = bofvar
	  chi2_i( i ) = chi2
	   bof_i( i ) = bof
	ndat_i( i ) = ndat
	   dof_i( i ) = dof
	   aic_i( i ) = aic
	   bic_i( i ) = bic
      end if

* KDH : SHOULD OMIT CALIBRATION PARAMETERS NOT USED
* KDH : ASSUMES ALL SCOPES OBSERVE AT ALL WAVELENGTHS :(

* calibration parameters (3 per scope per wavelength)
* assumes all scopes observed at all wavelengths
	dofcal = 3 * ns
	if( iti .le. 0 ) dofcal = 0
* wavelengths
      do iw = 1, nw
* gather
	ndat   =   ndat_w( iw )
	chi2   =    sum_w( iw )
	dofroa = dofroa_w( iw )
	bofvar = bofvar_w( iw )
* compute 
	dofdat = ndat
	doffit = dofcal + dofroa
	dof = dofdat - doffit
	bof = chi2 + bofvar
	aic = bof + doffit * 2.d0
	bic = bof + doffit * dlog( dofdat )
* stow
	  ndat_w( iw ) = ndat
	dofcal_w( iw ) = dofcal	
	bofvar_w( iw ) = bofvar
	   dof_w( iw ) = dof
	  chi2_w( iw ) = chi2
	   aic_w( iw ) = aic
	   bic_w( iw ) = bic
* next wavelength iw
      end do

* calibration parameters (3 per wavelength per scope)
* assumes all scopes observed at all wavelengths
	dofcal = 3 * nw
c	if( iter .le. 0 ) dofcal = 0
	if( iti .le. 0 ) dofcal = 0
* scopes
      do is = 1, ns
* gather
	ndat   =   ndat_s( is )
	chi2   =    sum_s( is )
	dofroa = dofroa_s( is )
	bofvar = bofvar_s( is )
* compute
	dofdat = ndat
	doffit = dofcal + dofroa
	dof = dofdat - doffit
	bof = chi2 + bofvar
	aic = bof + doffit * 2.d0
	bic = bof + doffit * dlog( dofdat )
* stow
	  ndat_s( is ) = ndat
	dofcal_s( is ) = dofcal
	   dof_s( is ) = dof
	  chi2_s( is ) = chi2
	bofvar_s( is ) = bofvar
	   aic_s( is ) = aic
	   bic_s( is ) = bic
* next is
      end do
	
* stow per iteration and scope
	i = iti
      if( i .gt. 0 ) then
      do is = 1, ns
	dofcal_is( i, is ) = dofcal_s( is )
	dofroa_is( i, is ) = dofroa_s( is )
	bofvar_is( i, is ) = bofvar_s( is )
	  chi2_is( i, is ) =   chi2_s( is )
	   dof_is( i, is ) =    dof_s( is )
	   aic_is( i, is ) =    aic_s( is )
	   bic_is( i, is ) =    bic_s( is )
      end do
* per iteration and wavelength
      do iw = 1, nw
	dofcal_iw( i, iw ) = dofcal_w( iw )
	dofroa_iw( i, iw ) = dofroa_w( iw )
	bofvar_iw( i, iw ) = bofvar_w( iw )
	  chi2_iw( i, iw ) =   chi2_w( iw )
	   dof_iw( i, iw ) =    dof_w( iw )
	   aic_iw( i, iw ) =    aic_w( iw )
	   bic_iw( i, iw ) =    bic_w( iw )
      end do

* per iteration, scope and wavelength
      do iw = 1, nw
      do is = 1, ns
	add_isw( i, is, iw ) = add_sw( is, iw )
	scl_isw( i, is, iw ) = scl_sw( is, iw )
	rms_isw( i, is, iw ) = rms_sw( is, iw )
	sigadd_isw( i, is, iw ) = sigadd_sw( is, iw )
	sigscl_isw( i, is, iw ) = sigscl_sw( is, iw )
	sigrms_isw( i, is, iw ) = sigrms_sw( is, iw )
      end do
      end do
      end if

* cpu time 
	call cpu( dcpu )
	cputime = cputime + dcpu
	cpudone = cpudone + dcpu
	write(*,*) 'dCPU', dcpu, ' CPU', cputime, ' CPU/it', cputime/its

* done and to go
	itsdone = its + last * ( ig - 1 )
	itstogo = ( last - its ) + last * ( ngrid - ig )
	cputogo = itstogo * cputime / its

* report done
	i = itsdone
	cpusec = cpudone
	cpumin = cpusec / 60.
	cpuhrs = cpumin / 60.
      if( cpusec .lt. 60. ) then
	write(*,*) 'Iterations done', i, ' CPU', cpusec, ' sec'
      else if( cpumin .lt. 60. ) then
	write(*,*) 'Iterations done', i, ' CPU', cpumin, ' min'
      else
	write(*,*) 'Iterations done', i, ' CPU', cpuhrs, ' hr'
      end if
* report togo
	i = itstogo
	cpusec = cputogo
	cpumin = cpusec / 60.
	cpuhrs = cpumin / 60.
      if( cpusec .lt. 60. ) then
	write(*,*) 'Iterations todo', i, ' CPU', cpusec, ' sec'
      else if( cpumin .lt. 60. ) then
	write(*,*) 'Iterations todo', i, ' CPU', cpumin, ' min'
      else
	write(*,*) 'Iterations todo', i, ' CPU', cpuhrs, ' hr'
      end if

	write(*,*) 'Grid', ig, ' HW', hw, ' Iter', iti, ' ---------', converged

* next iteration its
	if( converged ) exit
      end do

* cpu
      if( ngrid .gt. 1 ) then
	cpugrid = cpugrid + cputime
      else
	cpuiter = cpuiter + cputime
      end if

* report
      if( converged ) then
	write(*,*) 'Converged after iteration', iti, ' :)'
      else
	write(*,*) '** NOT CONVERGED after iteration', iti, ' :('
      end if

* stow grid results
      if( ngrid .gt. 1 ) then
* gather
	i = iti
	aic = aic_i( i )
	bic = bic_i( i )
	chi2 = chi2_i( i )
	dof = dof_i( i )
	bof = bof_i( i )
	bofvar = bofvar_i( i )	
	dofcal = dofcal_i( i )
	dofroa = dofroa_i( i )
	ndat = ndat_i( i )
* stow
	hw_g( ig ) = hw
	dt_g( ig ) = dt
	nits_g( ig ) = i
	ndat_g( ig ) = ndat
	aic_g( ig ) = aic
	bic_g( ig ) = bic
	chi2_g( ig ) = chi2
	dof_g( ig ) = dof
	bof_g( ig ) = bof
	bofvar_g( ig ) = bofvar
	dofcal_g( ig ) = dofcal
	dofroa_g( ig ) = dofroa
* calibrations
      do iw = 1, nw
      do is = 1, ns
	scl_gsw( ig, is, iw ) = scl_isw( i, is, iw )
	add_gsw( ig, is, iw ) = add_isw( i, is, iw )
	rms_gsw( ig, is, iw ) = rms_isw( i, is, iw )
	sigscl_gsw( ig, is, iw ) = sigscl_isw( i, is, iw )
	sigadd_gsw( ig, is, iw ) = sigadd_isw( i, is, iw )
	sigrms_gsw( ig, is, iw ) = sigrms_isw( i, is, iw )
      end do
      end do
* report
	write(*,*) 'Stowed grid', ig, ' of', ngrid, ' HW', hw
	write(*,*) 'Chi2', chi2, '-2lnL', bof
	write(*,*) ' AIC', aic, '  BIC', bic
      end if


* next grid point ig
      end do

	if( ngrid .le. 1 ) goto 11

* find best AIC and BIC
	iga = 1
	igb = 1
      do ig = 1, ngrid
	aic = aic_g( ig )
	bic = bic_g( ig )
	if( aic .lt. aic_g( iga ) ) iga = ig
	if( bic .lt. bic_g( igb ) ) igb = ig
      end do
	hwaic = hw_g( iga )
	hwbic = hw_g( igb )
	bestaic = aic_g( iga )
	bestbic = bic_g( igb )

* report
	write(*,*) 'Best AIC', bestaic, ' HW', hwaic
	write(*,*) 'Best BIC', bestbic, ' HW', hwbic

* parabolic fit to grid minimum
* aic min
      do i = 1, ngrid
	plotx( i ) = hw_g( i )
	ploty( i ) = aic_g( i )
      end do
	call fitmin( ngrid, plotx, ploty, i, hwaic, bestaic, slope, curv )
	sighwaic = -1.
	if( curv .gt. 0.d0 ) sighwaic = dsqrt( 2.d0 / curv )
* bic min
      do i = 1, ngrid
	plotx( i ) = hw_g( i )
	ploty( i ) = bic_g( i )
      end do
	call fitmin( ngrid, plotx, ploty, i, hwbic, bestbic, slope, curv )
	sighwbic = -1.
	if( curv .gt. 0.d0 ) sighwbic = dsqrt( 2.d0 / curv )
	write(*,*) 'Parabolic fit to 3 points:'
	write(*,*) 'Best AIC', bestaic, ' HW', hwaic, ' +/-', sighwaic
	write(*,*) 'Best BIC', bestbic, ' HW', hwbic, ' +/-', sighwbic

* adopt best BIC
      if( sighwbic .gt. 0. ) then
	hw = hwbic
	bic = bestbic
	write(*,*) 'Best BIC adopted'
* interpolate calibrations
	x = hwbic
	call lookup( ngrid, plotx, x, i1, i2, p )
      do iw = 1, nw
      do is = 1, ns
	scl = terp1( ngrid, scl_gsw( 1, is, iw ), i1, i2, p )
	add = terp1( ngrid, add_gsw( 1, is, iw ), i1, i2, p )
	rms = terp1( ngrid, rms_gsw( 1, is, iw ), i1, i2, p )
	scl_sw( is, iw ) = scl
	add_sw( is, iw ) = add
	rms_sw( is, iw ) = rms
	sigscl = terp1( ngrid, sigscl_gsw( 1, is, iw ), i1, i2, p )
	sigadd = terp1( ngrid, sigadd_gsw( 1, is, iw ), i1, i2, p )
	sigrms = terp1( ngrid, sigrms_gsw( 1, is, iw ), i1, i2, p )
	sigscl_sw( is, iw ) = sigscl
	sigadd_sw( is, iw ) = sigadd
	sigrms_sw( is, iw ) = sigrms
	write(*,*) 'Wave', iw, ' Scope', is
	write(*,*) ' scl', scl, ' add', add, ' rms', rms
	write(*,*) ' +/-', sigscl, ' +/-', sigadd, ' +/-', sigrms
      end do
      end do
	write(*,*) 'Calibrations interpolated.'
      end if

	goto 600


* ------------------------------------------
* grid scan plot
  650	write(*,*) '------------------------------'
	write(*,*) 'DREAM Grid Scan Plot'
* format
	ipanabc = 1
	ipanscl = 2
	ipanrms = 3
	ipanadd = 4
	nx = 2
	ny = 2
	npan = 4
* open
	call pgstart( nx, ny, nouse )
	call pgbgw
	csize = 1.5
	call pgsch( csize )

* KDH : WORKS ONLY FOR 1 FILTER -- GENERALISE
* filter
	iw = 1
* title
	title = 'DREAM iterations'
	call append_i( title, itg )
* x label
	xlabel = 'ROA'
	if( shape .eq. 'C' ) call append( xlabel, '1/Cosh' )
	if( shape .eq. 'L' ) call append( xlabel, 'Lorentzian' )
	if( shape .eq. 'G' ) call append( xlabel, 'Gaussian' )
	if( shape .eq. 'B' ) call append( xlabel, 'Boxcar' )
	call append( xlabel, 'halfwidth \gt' )
	write(*,*) xlabel(:len1(xlabel))

* x range ( same for all panels )
	nplot = ngrid
	if( ngrid .le. 1 ) nplot = nhw
      do i = 1, nplot
	plot1( i ) = hw_g( i )
      end do

	call mnmx( nplot, plot1, xmn, xmx )
	write(*,*) 'HW Grid', nplot, (plot1(i),i=1,nplot)
	write(*,*) ' x:', xmn, xmx
* log x
	logx = 0
	if( xmx .gt. 2. * xmn ) logx = 1
      if( logx .ne. 0 ) then
      do i = 1, nplot
	plot1(i) = alog10( plot1(i) ) 
      end do
	xmx = alog10( xmx )
	xmn = alog10( xmn )
	write(*,*) ' logx:', xmn, xmx
      end if
* expand x
	add = ( xmx - xmn ) * 0.05
	xmn = xmn - add
	xmx = xmx + add

* panels
      do ipan = 1, npan
	write(*,*) '---- Panel', ipan, ' of', npan, '-----'

* aic, bic, chi2
      if( ipan .eq. ipanabc ) then
	title = 'DREAM iterations'
	call append_i( title, itg )
	ylabel = 'Badness of Fit'
	write(*,*) ylabel(:len1(ylabel))
      do i = 1, nplot
	plot2( i ) = aic_g( i )
	plot3( i ) = bic_g( i )
	plot4( i ) = chi2_g( i )
* bof = -2lnL
	plot5( i ) = bof_g( i )
	plot6( i ) = dof_g( i )
* occam penalties
	plot7( i ) = aic_g( i ) - bof_g(i )
	plot8( i ) = bic_g( i ) - bof_g(i )
	plot9( i ) = bofvar_g( i )
      end do

* y limits
	call mnmx( nplot, plot2, amn, amx )
	call mnmx( nplot, plot3, bmn, bmx )
	call mnmx( nplot, plot4, cmn, cmx )
	call mnmx( nplot, plot5, emn, emx )
	call mnmx( nplot, plot6, dmn, dmx )
	ymx = max( amx, bmx, cmx, dmx, emx )
	ymn = min( amn, bmn, cmn, dmn, emn )
* report
	write(*,*) ' AIC ', amn, amx 
	write(*,*) ' BIC ', bmn, bmx 
	write(*,*) ' Chi2', cmn, cmx 
	write(*,*) ' DoF ', dmn, dmx 
	write(*,*) '-2lnL', emn, emx 
	write(*,*) ' y', ymn, ymx

* calibrations (per scope)
      else
* gather
* scl
      if( ipan .eq. ipanscl ) then
	ylabel = 'scl'
      do is = 1, ns
      do i = 1, nplot
	scl = scl_gsw( i, is, iw )
	sig = sigscl_gsw( i, is, iw )
	plot1_s( i, is ) = scl
	plot2_s( i, is ) = scl + sig
	plot3_s( i, is ) = scl - sig
      end do
      end do
* add
      else if( ipan .eq. ipanadd ) then
	ylabel = 'add'
      do is = 1, ns
      do i = 1, nplot
	add = add_gsw( i, is, iw )
	sig = sigadd_gsw( i, is, iw )
	plot1_s( i, is ) = add
	plot2_s( i, is ) = add + sig
	plot3_s( i, is ) = add - sig
      end do
      end do
* rms
      else if( ipan .eq. ipanrms ) then
	ylabel = 'rms'
      do is = 1, ns
      do i = 1, nplot
	rms = rms_gsw( i, is, iw )
	sig = sigrms_gsw( i, is, iw )
	plot1_s( i, is ) = rms
	plot2_s( i, is ) = rms + sig
	plot3_s( i, is ) = rms - sig
      end do
      end do
      end if
	write(*,*) ylabel(:len1(ylabel))

* y limits
      do i = 1, nplot
      do is = 1, ns	
	call mnmx( nplot, plot1_s(1,is), smn, smx )
c	write(*,*) 'Grid', i, smn, smx, ' Scope', is
      if( i .eq. 1 .and. is .eq. 1 ) then
	ymn = smn
	ymx = smx
      else
	ymx = max( ymx, smx )
	ymn = min( ymn, smn )
      end if
      end do
c	write(*,*) 'Grid', i, ymn, ymx
      end do
      end if

* log y
	logy = 0
      if( ipan .eq. ipanabc ) then
	if( ymn .gt. 0. .and. ymx .gt. 2. * ymn ) logy = 1
      if( logy .ne. 0 ) then
      do i = 1, nplot
	plot2( i ) = alog10( plot2( i ) ) 
	plot3( i ) = alog10( plot3( i ) ) 
	plot4( i ) = alog10( plot4( i ) ) 
	plot5( i ) = alog10( plot5( i ) ) 
	plot6( i ) = alog10( plot6( i ) ) 
	plot7( i ) = alog10( plot7( i ) ) 
	plot8( i ) = alog10( plot8( i ) ) 
	plot9( i ) = alog10( plot9( i ) ) 
      end do
	ymx = alog10( ymx )
	ymn = alog10( ymn )
	write(*,*) 'logy', ymn, ymx
      end if
      end if
* expand
	add = ( ymx - ymn ) * 0.05
	ymn = ymn - add
	ymx = ymx + add
* window
	logxy = 10 * logx + 20 * logy
	write(*,*) 'pgenv(', xmn, xmx, ymn, ymx, logxy, ' )'
	call pgenv( xmn, xmx, ymn, ymx, 0, logxy )
* grid
	call pgsls( 4 )
c	call pgbox( 'GL', 0.0, 0, 'GL', 0.0, 0 )
	if( logx .ne. 0 ) call pgbox( 'GL', 1., 0, ' ', 0., 0 )
	if( logy .ne. 0 ) call pgbox( ' ', 1., 0, 'GL', 0., 0 )

	call pgsls( 1 )
* labels
	call pglabel( xlabel, ylabel, title )

* data and filter
	label = ' '
      if( obsfile .ne. ' ' .and. obsfile .ne. 'FAKE') then
	call append( label, obsfile )
      end if
	call append( label, 'filter:' )
	call append( label, name_w( iw ) )
	call pgmtxt( 'T', 0.5, 0., 0., label )
* scope(s)
      if( ipan .ne. ipanabc ) then
	label = ' '
	call append_i( label, ns )
	call pgmtxt( 'T', 0.5, 1., 1., label )
	label = ' '
	call append( label, ' scopes:' )
	call pgmtxt( 'T', 0.5, 1., 0., label )
      end if

* verticals at best aic and bic
      do iab = 1, 2
      if( iab .eq. 1 ) then
	call pgsci( kred )
	cen = hwaic
	sig = sighwaic
      else
	call pgsci( kblue )
	cen = hwbic
	sig = sighwbic
      end if
	kk = 0
	if( sig .gt. 0. ) kk = 1
      do k = -kk, kk
	x = cen + k * sig
	if( logx .ne. 0 ) x = alog10( x )
	call pgsls( 1 + iabs( k ) )
	call pgmove( x, ymn )
	call pgdraw( x, ymx )
	call pgsls( 1 )
      end do
      end do

	ibest = istar
	isymb = icirc
      if( ipan .eq. ipanabc ) then
* aic -----------------------
	call pgsci( kred )
	call pgpoint( nplot, plot1, plot2, isymb )
	call pgline( nplot, plot1, plot2 )

	call pgsls( 2 )
	call pgpoint( nplot, plot1, plot7, isymb )
	call pgline( nplot, plot1, plot7 )
	call pgsls( 1 )


	call pgmtxt( 'T', -3.0, 1., 0., ' AIC' )
	write(*,*) 'nplot', nplot
	write(*,*) ' HW ', (plot1(i),i=1,nplot)
	write(*,*) ' AIC', (plot2(i),i=1,nplot)
* best
	write(*,*) 'Best AIC', bestbic, ' HW', real(hwaic)
     &		, '+/-', real(sighwaic)
	x = hwaic
	y = bestaic
	if( logx .ne. 0 ) x = alog10( x )
	if( logy .ne. 0 ) y = alog10( y )
	write(*,*) x, y
	call pgpoint( 1, x, y, ibest )
* bic -----------------------
	call pgsci( kblue )
	call pgpoint( nplot, plot1, plot3, isymb )
	call pgline( nplot, plot1, plot3 )

	call pgsls( 2 )
	call pgpoint( nplot, plot1, plot8, isymb )
	call pgline( nplot, plot1, plot8 )
	call pgsls( 1 )

	call pgmtxt( 'T', -1.5, 1., 0., ' BIC' )
	write(*,*) ' BIC', (plot3(i),i=1,nplot)
* best
	write(*,*) 'Best BIC', bestbic, ' HW', real(hwbic)
     &		, '+/-', real(sighwbic)
	x = hwbic
	y = bestbic
	if( logx .ne. 0 ) x = alog10( x )
	if( logy .ne. 0 ) y = alog10( y )
	write(*,*) x, y
	call pgpoint( 1, x, y, ibest )
* BoF -----------------------
	call pgsci( kblack )
	call pgpoint( nplot, plot1, plot5, isymb )
	call pgline ( nplot, plot1, plot5 )
	call pgmtxt( 'T', -4.5, 1.0, 0., ' -2lnL' )

	call pgsls( 2 )
	call pgpoint( nplot, plot1, plot9, isymb )
	call pgline( nplot, plot1, plot9 )
	call pgsls( 1 )

* chi2 ----------------------
	call pgsci( kpurple )
	call pgpoint( nplot, plot1, plot4, isymb )
	call pgline ( nplot, plot1, plot4 )
	call pgmtxt( 'T', -7.5, 1., 0., ' \gx\u2\d' )
	write(*,*) 'Chi2', (plot4(i),i=1,nplot)
* DoF -----------------------
	call pgsci( kgreen )
	call pgpoint( nplot, plot1, plot6, isymb )
	call pgline ( nplot, plot1, plot6 )
	call pgmtxt( 'T', -9.0, 1.0, 0., ' DoF' )
	write(*,*) ' DoF', (plot5(i),i=1,nplot)

* calibration per scope
      else

* horizontal at scl = 1.
	call pgsci( kblack )
	call pgsls( 4 )
      if( ipan .eq. ipanscl ) then
	y = 1.
	if( logy .ne. 0 ) y = 0.
	call pgmove( xmn, y )
	call pgdraw( xmx, y )
* horizontal at add = 0
      else if( ipan .eq. ipanadd ) then
	y = 0.
	call pgmove( xmn, y )
	call pgdraw( xmx, y )
      end if
	call pgsls( 1 )


      do is = 1, ns
	call pgsci( kolor_s( is ) )
	call pgpoint( nplot, plot1, plot1_s(1,is), isymb )
	call pgline ( nplot, plot1, plot1_s(1,is) )
	call pgsls( 2 )
	call pgline ( nplot, plot1, plot2_s(1,is) )
	call pgline ( nplot, plot1, plot3_s(1,is) )
	call pgsls( 1 )
	name = name_s( is )
	y = -1.5 * is
	call pgmtxt( 'T', y, 1., 0., ' ' // name )

* star 
	x = hw
	if( logx .ne. 0 ) x = alog10( x )
      if( ipan .eq. ipanscl ) then
	y = scl_sw( is, iw )
      else if( ipan .eq. ipanadd ) then
	y = add_sw( is, iw )
      else if( ipan .eq. ipanrms ) then
	y = rms_sw( is, iw )
      end if
	if( logy .ne. 0 ) y = alog10( y )
	write(*,*) 'Wave', iw, ' Scope', is, x, y
	call pgsch( 1.5 * csize )
	call pgpoint( 1, x, y, istar )
	call pgsch( csize )

* next scope
      end do
	call pgsci( kblack )

      end if

	call pgsci( kblack )

* next panel
      end do

* close
	call pgstop
* again
	if( pgagain() ) goto 650
	goto 600	

* ------------------------------------------
* lightcurve plot

* format
  100 if( nw .eq. 1 ) then
	nx = 1
	ny = 3
	npan = nx * ny
	irowdat = 1
	irowcal = 2
	irowres = 3
      else
	nx = 2
	ny = nw
	icoldat = 0
	icolcal = 1
	icolres = 2
      end if
	npan = nx * ny
* open
	call pgstart( nx, ny, nouse )
	call pgbgw
	csize = 1.5
	call pgsch( csize )

* evaluate current model and residuals

      do io = 1, no
	 iw =  iw_o( io )
	 is =  is_o( io )
	dat = dat_o( io )
	  t =   t_o( io )

	nt = nt_w( iw )
	call lookup( nt, t_tw( 1, iw ), t, i1, i2, p )
	x = terp1(   nt, x_tw( 1, iw ),    i1, i2, p )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	fit = add + scl * x
* stow
	fit_o( io ) = fit
	dif_o( io ) = dat - fit
      end do

* calibrated data

      do io = 1, no
	dat = dat_o( io )
	sig = sig_o( io )
	dif = dif_o( io )
	iw = iw_o( io )
	is = is_o( io )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
* old error bar
	sigold = sig
* new error bar
	sig = sqrt( sig * sig + rms * rms )

	cal = ( dat - add ) / scl
	difcal = dif / scl
	sigcal = sig / scl
	rmscal = rms / scl
	sigoldcal = sigold / scl
* stow
	   cal_o( io ) = cal
	difcal_o( io ) = difcal
	sigcal_o( io ) = sigcal
	rmscal_o( io ) = rmscal
	sigoldcal_o( io ) = sigoldcal
      end do

* panels
      do ipan = 1, npan
	irow = 1 + ( ipan - 1 ) / nx
	icol = ipan - nx * ( irow - 1 )
	write(*,*) 'panel', ipan, ' col', icol, ' row', irow
* what to plot
	ipandat = 0
	ipancal = 0
	ipanres = 0
      if( nw .eq. 1 ) then
	iw = 1
	if( irow .eq. irowcal ) ipancal = ipan
	if( irow .eq. irowdat ) ipandat = ipan
	if( irow .eq. irowres ) ipanres = ipan
      else
	iw = irow
	if( icol .eq. icolcal ) ipancal = ipan
	if( icol .eq. icoldat ) ipandat = ipan
	if( icol .eq. icolres ) ipanres = ipan
      end if

* gather for wavelength iw data from all scopes
      do is = 1, ns
	n = 0
      do io = 1, no
      if( iw_o( io ) .eq. iw ) then
      if( is_o( io ) .eq. is ) then
      if( n .lt. maxt ) then
	n = n + 1
	  t_ts( n, is ) =   t_o( io )
	dat_ts( n, is ) = dat_o( io )
	sig_ts( n, is ) = sig_o( io )

	   cal_ts( n, is ) =    cal_o( io )
	difcal_ts( n, is ) = difcal_o( io )
	sigcal_ts( n, is ) = sigcal_o( io )
	rmscal_ts( n, is ) = rmscal_o( io )
	sigoldcal_ts( n, is ) = sigoldcal_o( io )

	fit_ts( n, is ) = fit_o( io )
	dif_ts( n, is ) = dif_o( io )
	chi_ts( n, is ) = chi_o( io )
      end if
      end if
      end if
      end do
	nt_s( is ) = n
* next scope
      end do

* gather calibrations for this wavelength
      do is = 1, ns
	add_s( is ) = add_sw( is, iw )
	scl_s( is ) = scl_sw( is, iw )
	rms_s( is ) = rms_sw( is, iw )
      end do

* gather model lightcurve for this wavelength
	nt = nt_w( iw )
	write(*,*) ' iw', iw, ' nt_w', nt
      do it = 1, nt
	   t_t( it ) =    t_tw( it, iw )
	   x_t( it ) =    x_tw( it, iw )
	sigx_t( it ) = sigx_tw( it, iw )
      end do

* x label
	xlabel = 'time'
	title = 'DREAM test on FAKE data'
      if( obsfile .ne. ' ' .and. obsfile .ne. 'FAKE') then
	title = 'DREAM fit to'
	call append( title, obsfile )
	call append( title, 'data' )
      if( jd0 .ne. 0 ) then
	xlabel = 'HJD -'
	call append_I( xlabel, jd0 )
      end if
      end if
* title
	call append( title, name_w( iw ) )
	lam = nint( w_w( iw ) )
      if( lam .gt. 0 ) then
	call append_i( title, lam )
	call append( title, '\A' )
      end if


* -------------------------------------------
      if( ipan .eq. ipandat ) then
* x range
	call mnmx( nt, t_t, xmn, xmx )
	write(*,*) 'Times', nt, xmn, xmx 

* y range
c	ymn = 0.
c	ymx = 0.
      do is = 1, ns
	n = nt_s( is )
	call mnmx( n, dat_ts(1,is), dmn, dmx )
c	call mnmx( n, dif_ts(1,is), emn, emx )
c	ymn = min( ymn, dmn, emn )
c	ymx = max( ymx, dmx, emx )
      if( is .eq. 1 ) then
	ymn = dmn
	ymx = dmx
      else
	ymn = min( ymn, dmn )
	ymx = max( ymx, dmx )
      end if
      end do
* expand x
	add = ( xmx - xmn ) * 0.1
	xmn = xmn - add
	xmx = xmx + add
* expand y
	add = ( ymx - ymn ) * 0.2
	ymx = ymx + add
	ymn = ymn - add
* window
	call pgenv( xmn, xmx, ymn, ymx, 0, 0 )
	call pgmove( xmn, 0. )
	call pgdraw( xmx, 0. )
* grid
	call pgsls( 4 )
	call pgbox( 'G', 0.0, 0, 'G', 0.0, 0 )
	call pgsls( 1 )
* labels
	ylabel = 'flux data'
	call pglabel( xlabel, ylabel, title )

* sort scopes by mean variance
      do is = 1, ns
	n = nt_s( is )
      do i = 1, n
	sig = sig_ts( it, is )
	plot1( i ) = sig * sig
      end do
	call avgrms( n, plot1, avg, rms )
	plot2( is ) = avg
      end do
	call quicksort( ns, plot2, i_s )

* scopes with larger errors first, then smaller ones
      do iis = ns, 1, -1
	is = i_s( iis )

	isymb = idot
	call pgsci( kolor_s( is ) )

* KDH : OMIT MODELS ON RAW FLUX PLOT
      IF( .FALSE. ) THEN
	rms = rms_s( is )
	add = add_s( is )
	scl = scl_s( is )
      do it = 1, nt
	x = x_t( it )
	fit = add + scl * x
	fit_t( it ) = fit
      end do
* error envelope
      do k = -1, 1
      do it = 1, nt
	plot1( it ) = fit_t( it ) + k * rms
      end do
	call pgline( nt, t_t, plot1 )
      end do
      END IF

* data
	n = nt_s( is )
	call pgerrb( 6, n, t_ts(1,is), dat_ts(1,is), sig_ts(1,is), term )
	call pgpoint(   n, t_ts(1,is), dat_ts(1,is), isymb )
* residuals
      IF( .FALSE. ) THEN
	call pgerrb( 6, n, t_ts(1,is), dif_ts(1,is), sig_ts(1,is), term )
	call pgpoint(   n, t_ts(1,is), dif_ts(1,is), isymb )
      END IF
* next scope is
      end do
	call pgsci( kblack )

* -------------------------------------------
* calibrated data
      else if( ipan .eq. ipancal .or. ipan .eq. ipanres ) then
* x range
	call mnmx( nt, t_t, xmn, xmx )
* y range
      do is = 1, ns
	n = nt_s( is )
      if( ipan .eq. ipancal ) then
	call mnmx( n,    cal_ts(1,is), dmn, dmx )
      if( is .eq. 1 ) then
	ymn = dmn
	ymx = dmx
      else
	ymn = min( ymn, dmn )
	ymx = max( ymx, dmx )
      end if
* residuals
      else if( ipan .eq. ipanres ) then
	call mnmx( n, difcal_ts(1,is), emn, emx )
      if( is .eq. 1 ) then
	ymn = min( 0., emn )
	ymx = max( 0., emx )
      else
	ymn = min( ymn, emn )
	ymx = max( ymx, emx )
      end if
      end if
* next scope is
      end do
* expand x
	add = ( xmx - xmn ) * 0.1
	xmn = xmn - add
	xmx = xmx + add
* expand y
	add = ( ymx - ymn ) * 0.2
	ymx = ymx + add
	ymn = ymn - add
* window
	call pgenv( xmn, xmx, ymn, ymx, 0, 0 )

* grid
	call pgsls( 4 )
	call pgbox( 'G', 0.0, 0, 'G', 0.0, 0 )
	call pgsls( 1 )
* horizontal
	call pgmove( xmn, 0. )
	call pgdraw( xmx, 0. )

* labels
	ylabel = 'inter-calibrated flux data'
	if( ipan .eq. ipanres ) ylabel = 'residuals'
	call pglabel( xlabel, ylabel, title )


* sort scopes by mean variance
      do is = 1, ns
	n = nt_s( is )
      do i = 1, n
	sig = sigcal_ts( it, is )
	plot1( i ) = sig * sig
      end do
	call avgrms( n, plot1, avg, rms )
	plot2( is ) = avg
      end do
	call quicksort( ns, plot2, i_s )

* data first then model, or vice versa
	mfit = 1
	mdat = 2
      do m=1,2
      if( m .eq. mdat ) then

* scopes with larger errors first, then smaller
      do iis = ns, 1, -1
	is = i_s( iis )
	n = nt_s( is )

	isymb = idot
	call pgsci( kolor_s( is ) )

* original calibrated error bars
      if( ipan .eq. ipancal .or. ipan .eq. ipanres ) then
      do i = 1, n
	sigold = sigoldcal_ts( i, is )
	plot1( i ) = sigold
      end do
      end if

* calibrated lightcurve data
      if( ipan .eq. ipancal ) then
* original calibrated error bars
	call pgerrb( 6, n, t_ts(1,is), cal_ts(1,is), plot1, term )
* revised calibrated error bars
	call pgerrb( 6, n, t_ts(1,is), cal_ts(1,is), sigcal_ts(1,is), term )
	call pgpoint(   n, t_ts(1,is), cal_ts(1,is), isymb )

* residuals of calibrated lightcurve data
      else if( ipan .eq. ipanres ) then
* original calibrated error bars
	call pgerrb( 6, n, t_ts(1,is), difcal_ts(1,is), plot1, term )
* revised calibrated error bars
	call pgerrb( 6, n, t_ts(1,is), difcal_ts(1,is), sigcal_ts(1,is), term )
	call pgpoint(   n, t_ts(1,is), difcal_ts(1,is), isymb )
      end if
      end do	

      else if( m .eq. mfit ) then

* model lightcurve and error envelope
      if( ipan .eq. ipancal ) then
	call pgsci( kblack )
      do k = -1, 1
      do it = 1, nt
	plot1( it ) = x_t( it ) + k * sigx_t( it )
      end do
	call pgsci( kblack )
	if( k .ne. 0 ) call pgsci( kdkgrey )
	call pgline( nt, t_t, plot1 )
      end do
      end if

* next m
      end if
      end do

	call pgsci( kblack )

* last panel
      end if

* ROA window function
      if( ipan .eq. ipancal .or. ipan .eq. ipanres ) then
      if( hw .gt. 0. .and. iter .gt. 0 ) then
	p = 0.5
	t1 = xmn * ( 1. - p ) + p * xmx
	call roa( shape, hw, 1, t1, 0., 1., 1, t, dat1, sig1 )
      do it = 1, nt
	t = t_t( it )
	call roa( shape, hw, 1, t1, dat1, sig1, 1, t, avg, sig )
	w = 0.
	if( sig .gt. 0. ) w = sig1 / sig
	win_t( it ) = w * w
      end do
	y = 3.
	call pgsci( kblue )
	call pgswin( xmn, xmx, 0., y )
	call pgbin( nt, t_t, win_t, .true. )
c	call pgline( nt, t_t, win_t )
	call pgswin( xmn, xmx, ymn, ymx )

	label = 'ROA ' // shape(:len1(shape))
	call append_n( label, hw, 3 )
	y = -0.5
	x = 0.03
	call pgmtxt( 'B', y, x, 0., label )
	call pgsci( kblack )
      end if
      end if

* scopes
	label = ' scope:'
      if( ns .gt. 1 ) then
	label = ' '
	call append_i( label, ns )
	call pgmtxt( 'T', 0.5, 1., 1., label )
	label = ' scopes:'
      end if
	call pgsci( kblack )
	call pgmtxt('T', 0.5, 1.0, 0., label ) 

	dy = -1.5

	call mnmx( nt, t_t, tmn, tmx )
	call mnmx( no, t_o, tomn, tomx )
      do is = 1, ns
	name = name_s( is )
	call pgsci( kolor_s( is ) )
	y = is * dy
	call pgmtxt( 'T', y, 1., 0., ' ' // name )
	call append( title, name_s( is ) )

* systematic error bars
	rms = rms_sw( is, iw )
	scl = scl_sw( is, iw )
	if( ipan .eq. ipancal ) rms = rms * scl
	p = real( is ) / max( 1, ns + 2 )
	x = tmx * ( 1. - p ) + p * xmx
	y = 0.
	call pgerrb( 6, 1, x, y, rms, term )
* next scope
      end do
	call pgsci( kblack )

* recall
	i = iter
	dofcal = dofcal_iw( i, iw )
	dofroa = dofroa_iw( i, iw )
	bofvar = bofvar_iw( i, iw )
	  chi2 =   chi2_iw( i, iw )
	   dof =    dof_iw( i, iw )
	   aic =    aic_iw( i, iw )
	   bic =    bic_iw( i, iw )
	ndat = ndat_w( iw )
* compute
	bof = chi2 + bofvar
	doffit = dofcal + dofroa
	chi2red = chi2 / max( 1., dof )

* report (warning)
	write(*,*) 'Ndat', ndat, ' Npar', doffit, 'Ndof', dof
      if( ndat .ne. dof + doffit ) then
	write(*,*) '** WARNING : Ndat - Npar - Ndof', n - doffit - dof
      end if
	write(*,*) 'Fit parameters: cal', dofcal, ' ROA', dofroa
	write(*,*) 'DoF =', n, ' -', doffit, ' =', dof
	write(*,*) 'Chi2', chi2, ' bofvar', bofvar
	write(*,*) '-2lnL', bof, ' AIC', aic, ' BIC', bic

* Chi^2/dof label	
      if( ipan .eq. ipancal .or. ipan .eq. ipanres ) then
	label = '\gx\u2\d='
	call append_f( label, real(chi2), '(f15.1)' )
	call append( label, '/' )
      if( doffit .le. 0. ) then
	call append_i( label, ndat )
      else
	call append( label, '(' )
	call append_i( label, ndat )
      if( dofcal .gt. 0. ) then
	call append( label, '-' )
	call append_i( label, nint(dofcal) )
      end if
      if( dofroa .gt. 0. ) then
	call append( label, '-' )
	call append_f( label, real(dofroa), '(F15.1)' )
      end if
	call append( label, ')' )
      end if
	call append( label, '=' )
	call append_n( label, chi2red, 3 ) 
	call nospace( label )
* bofvar
      if( bofvar .ne. 0. ) then
	call append( label, '\gSln\gs\u2\d' )
	call append_f( label, real(bofvar), '(f15.1)' )
      end if	

	call pgmtxt( 'T', 0.5, 0.0, 0.0, label )

* chi2, AIC, BIC
	label = ' '
	call append( label, '-2lnL' )
	call append_f( label, real(chi2 + bofvar), '(f15.1)' )
	call append( label, 'AIC' )
	call append_f( label, real(aic), '(f15.1)' )
	call append( label, 'BIC' )
	call append_f( label, real(bic), '(f15.1)' )
	call pgmtxt( 'T', -1.5, 0.05, 0., label )
      end if

* next panel
      end do
* close
	call pgstop
* again
	if( pgagain() ) goto 100
	goto 11

* convergence plot ----------------------------------

  200	n = iter
      if( n .lt. 1 ) then
	write(*,*) '** No Iterations yet.'
	goto 11
      end if
* format
	nx = 1
	ny = 2
	npan = 2
	ipanstep = 1
	ipanabc = 2

* open
	call pgstart( nx, ny, nouse )
	call pgbgw
	csize = 1.5
	call pgsch( csize )
* panels
      do ipan = 1, npan

* log x option
	logx = 0
	if( reply .eq. 'X' ) logx = 1

* load data for plot
      if( ipan .eq. ipanabc ) then
	logy = 1
	floor = 1.d-5
      do i = 1, n
	x = i
	a = aic_i( i )
	b = bic_i( i )
	c = chi2_i( i )
	d = dof_i( i )
	e = bof_i( i )
	if( logx .ne. 0 ) x = alog10( x )
	if( logy .ne. 0 ) a = alog10( max( floor, a ) )
	if( logy .ne. 0 ) b = alog10( max( floor, b ) )
	if( logy .ne. 0 ) c = alog10( max( floor, c ) )
	if( logy .ne. 0 ) d = alog10( max( floor, d ) ) 
	if( logy .ne. 0 ) e = alog10( max( floor, e ) ) 
	plot1( i ) = x
	plot2( i ) = a
	plot3( i ) = b
	plot4( i ) = c
	plot5( i ) = d
	plot6( i ) = e
      end do

      else if( ipan .eq. ipanstep ) then
	logy = 1
* load data for plot
      do i = 1, n
	chi = 0.
      do j = 1, nstep
	chi = max( chi, chi_ic(i,j) )
      end do
	x = i
	y = chi
	if( logx .ne. 0 ) x = alog10( x )
	if( logy .ne. 0 ) y = alog10( y )
	plot1( i ) = x
	plot2( i ) = y
* safe
	y = safe_i( i )
	if( logy .ne. 0 ) y = alog10( y )
	plot4( i ) = y
* tiny
	y = tiny_i( i )
	if( logy .ne. 0 ) y = alog10( y )
	plot5( i ) = y
* chi2
	ndof = dof_i( i )
	y = chi2_i( i ) / max( 1, ndof )
	if( logy .ne. 0 ) y = alog10( y )
	plot6( i ) = y
* dChi2
	floor = tiny / 100.
	y = chi2_i( i )
	if( i .gt. 1 ) y = chi2_i( i ) - chi2_i( i-1 )
	if( logy .ne. 0 ) y = alog10( max( floor, abs(y) ) )
	plot7( i ) = y 
* dAIC
	y = aic_i( i )
	if( i .gt. 1 ) y = aic_i( i ) - aic_i( i-1 )
	if( logy .ne. 0 ) y = alog10( max( floor, abs(y) ) )
	plot8( i ) = y 
* dBIC
	y = bic_i( i )
	if( i .gt. 1 ) y = bic_i( i ) - bic_i( i-1 )
	if( logy .ne. 0 ) y = alog10( max( floor, abs(y) ) )
	plot9( i ) = y 

* next iteration i
      end do

* end if ipanstep
      end if


* x range
	xmn = 0.5
	xmx = n + 0.5
      if( logx .ne. 0 ) then      
	xmn = alog10( xmn )
	xmx = alog10( xmx )
      end if

	add = ( xmx - xmn ) * 0.05
	xmn = xmn - add
	xmx = xmx + add

* y range
* badness of fit
      if( ipan .eq. ipanabc ) then
	call mnmx( n, plot2, amn, amx )
	call mnmx( n, plot3, bmn, bmx )
	call mnmx( n, plot4, cmn, cmx )
	call mnmx( n, plot5, dmn, dmx )
	call mnmx( n, plot6, emn, emx )
	ymn = min( amn, bmn, cmn, dmn, emn )
	ymx = max( amx, bmx, cmx, dmx, emx ) 

* parameter steps
      else if( ipan .eq. ipanstep ) then
	call mnmx( n, plot2, ymn, ymx )
c	call mnmx( n, plot4, ymn, ymx )
      do j = 1, nstep
      do i = 1, n
	y = chi_ic( i, j )
	if( logy .ne. 0 ) y = alog10( y )
	plot3( i ) = y
      end do
	call mnmx( n, plot3, cmn, cmx )
	ymn = min( ymn, cmn )
	ymx = max( ymx, cmx )
      end do
	y = tiny
	if( logy .ne. 0 ) y = alog10( y )
	ymn = max( ymn, y )
* end if ipanstep
      end if

* expand y range
	add = ( ymx - ymn ) * 0.05
	ymn = ymn - add
	ymx = ymx + add

* window
	logxy = 10 * logx + 20 * logy
	if( logy .ne. 0 .and. ymx - ymn .gt. 10 ) logxy = 10 * logx
	call pgenv( xmn, xmx, ymn, ymx, 0, logxy )
* grid
	call pgsls( 4 )
	if( logx .ne. 0 ) call pgbox( 'G', 1., 0, ' ', 0., 0 )
	if( logy .ne. 0 ) call pgbox( ' ', 1., 0, 'G', 0., 0 )
	call pgsls( 1 )

* horizontal
	y = tiny
	if( logy .ne. 0 ) y = alog10( y )
	call pgsls( 4 )
	call pgsci( kyellow )
	call pgmove( xmn, y )
	call pgdraw( xmx, y )
	call pgsls( 1 )
	call pgsci( kblack )

* labels
	xlabel = 'iterations'
	call append_i( xlabel, n )
      if( ipan .eq. ipanabc ) then
c	ylabel = 'AIC, BIC, \gx\u2\d'
	ylabel = 'Badness of Fit'
      else
	ylabel = 'step / sigma'
      end if
      if( logxy .ne. 10 * logx + 20 * logy ) then
	ylabel = 'log[ ' // ylabel(:len1(ylabel)) //' ]'
      end if
	title = 'DREAM Convergence Test'
	call pglabel( xlabel, ylabel, title )

 
* historams
	n10 = min( n, 10 )
	isymb = icirc
      if( ipan .eq. ipanabc ) then
* aic
	call pgsci( kred )
	call pgpoint( n10, plot1, plot2, isymb )
	call pgbin( n, plot1, plot2, .true. )
c	y = plot2( n )
c	call pgtext( x, y, 'AIC' )
	call pgmtxt( 'T', -3.0, 1., 0., ' AIC' )
* bic
	call pgsci( kblue )
	call pgpoint( n10, plot1, plot3, isymb )
	call pgbin( n, plot1, plot3, .true. )
c	y = plot3( n )
c	call pgtext( x, y, 'BIC' )
	call pgmtxt( 'T', -1.5, 1., 0., ' BIC' )
* bof
	call pgsci( kblack )
	call pgpoint( n10, plot1, plot6, isymb )
	call pgbin( n, plot1, plot6, .true. )
c	y = plot6( n )
c	call pgtext( x, y, '-2lnL' )
	call pgmtxt( 'T', -4.5, 1., 0., ' -2lnL' )
* chi2
	call pgsci( kpurple )
	call pgpoint( n10, plot1, plot4, isymb )
	call pgbin( n, plot1, plot4, .true. )
	x = xmx + 0.01 * ( xmx - xmn )
c	y = plot4( n )
c	call pgtext( x, y, '\gx\u2\d' )
	call pgmtxt( 'T', -7.5, 1., 0., ' \gx\u2\d' )

* dof
	call pgsci( kgreen )
	call pgpoint( n10, plot1, plot5, isymb )
	call pgbin( n, plot1, plot5, .true. )
c	y = plot5( n )
c	call pgtext( x, y, '\gSln\gs\u2\d' )
c	call pgtext( x, y, 'dof' )
	call pgmtxt( 'T', -9.0, 1., 0., ' dof' )

* step / sigma --------------
      else if( ipan .eq. ipanstep ) then
	call pgsci( kblue )
	call pgpoint( n10, plot1, plot2, isymb )
	call pgbin( n, plot1, plot2, .true. )

* safe
	call pgsci( korange )
	call pgpoint( n10, plot1, plot4, isymb )
	call pgsls( 2 )
	call pgbin( n, plot1, plot4, .true. )
	call pgsls( 1 )
	y = plot4( 1 )
	x = xmn
	call pgtext( x, y, ' safe' )
* tiny
	call pgsci( korange )
	call pgpoint( n10, plot1, plot5, isymb )
	call pgsls( 2 )
	call pgbin( n, plot1, plot5, .true. )
	call pgsls( 1 )
	y = plot5( 1 )
	x = xmn
	if( y .gt. ymn ) call pgtext( x, y, ' tiny' )
* chi2
	call pgsci( kpurple )
	call pgpoint( n10, plot1, plot6, isymb )
	call pgbin( n, plot1, plot6, .true. )
	y = plot6( n )
	call pgtext( xmx, y, ' \gx\u2\d/dof' )



	dy = -1.5
	y = - dy
* dchi2
	call pgsci( kpurple )
	call pgpoint( n10, plot1, plot7, isymb )
	call pgbin( n, plot1, plot7, .true. )
c	y = plot7( n )
c	call pgtext( xmx, y, ' \gD\gx\u2\d' )
	y = y + dy
	call pgmtxt( 'B', y, 1., 0., ' \gD\gx\u2\d' )

* dAIC
	call pgsci( kred )
	call pgpoint( n10, plot1, plot8, isymb )
	call pgbin( n, plot1, plot8, .true. )
d	y = plot8( n )
d	call pgtext( xmx, y, ' \gDAIC' )
	y = y + dy
	call pgmtxt( 'B', y, 1., 0., ' \gDAIC' )

* dBIC
	call pgsci( kblue )
	call pgpoint( n10, plot1, plot9, isymb )
	call pgbin( n, plot1, plot9, .true. )
c	y = plot9( n )
c	call pgtext( xmx, y, ' \gDBIC' )
	y = y + dy
	call pgmtxt( 'B', y, 1., 0., ' \gDBIC' )


* steps
      do j = 1, nstep
        name = name_c( j )
      do i = 1, n
        y = chi_ic( i, j )
        if( logy .ne. 0 ) y = alog10( y )    
        plot3( i ) = y
      end do
        call pgsci( j )
        call pgpoint( n10, plot1, plot3, isymb )
        call pgbin( n, plot1, plot3, .true. )
* legend
        dy = 1.5
        y = - dy * j
        x = 0.95
        call pgmtxt( 'T', y, x, 1., '\gD'//name )
      end do

* end if ipanstep
      end if
	call pgsci( kblack )

	i = iter
c	ndof = nint( dof_i( i ) )
c	npar = nint( dofcal_i( i ) )
c	chi2red = chi2_i( i ) / max( 1, ndof ) 
c	label = '\gx\u2\d/('
c	call append_i( label, ndof )
c	call append( label, '-' )
c	call append_i( label, npar )
c	call append( label, ')=' )
c	call append_n( label, chi2red, 3 ) 
c	call nospace( label )

* reduced Chi^2/dof
	ndat = ndat_i( i )
	chi2 = chi2_i( i )
	dofcal = dofcal_i( i )
	dofroa = dofroa_i( i )
	doffit = dofcal + dofroa
	dof = dof_i( i )
c	dof = ndof - dofcal - dofroa
	chi2 = chi2_i( i )
	chi2red = chi2 / max( 1.d0, dof )

* Chi^2/dof label	
	label = '\gx\u2\d='
	call append_f( label, real(chi2), '(f15.1)' )
	call append( label, '/' )
      if( doffit .le. 0. ) then
	call append_i( label, nint( doffit ) )
      else
	call append( label, '(' )
	call append_i( label, ndat )
      if( dofcal .gt. 0. ) then
	call append( label, '-' )
	call append_i( label, nint(dofcal) )
      end if
      if( dofroa .gt. 0. ) then
	call append( label, '-' )
	call append_f( label, real(dofroa), '(F15.1)' )
      end if
	call append( label, ')' )
      end if
	call append( label, '=' )
	call append_n( label, chi2red, 3 ) 
	call nospace( label )


	call pgmtxt( 'T', 0.5, 0.03, 0., label )

* next panel ipan
      end do
* close
	call pgstop
* again
	if( pgagain() ) goto 200
	goto 11

* load data from files
  500	write(*,*)
	write(*,*) 'Default extension is .lis'
	call spawn( 'ls -l *.lis' )
	write(*,*) 'Input file holding filter scope datafile'
	call inq( 'File', obsfile, obsfile )
	call word1( reply, l1, l2 )
	reply = reply( l1 : l2 )

	nf = maxf
	call loadlist( '@'//obsfile, nf, file_f )
	write(*,*)
	write(*,*) 'Observation Files', nf
	if( nf .le. 0 ) goto 11

* reset
	no = 0
	nw = 0
	ns = 0
	iter = 0
	
      do i = 1, nf
	file = file_f( i )
	write(*,*) file(:len1(file))
* filter
	call word1( file, l1, l2 )
	name = file(l1:l2)
	file = file(l2+1:)
	iw = 0
      do j = 1, nw
	if( name .eq. name_w( j ) ) iw = j
	if( iw .ne. 0 ) exit
      end do
      if( iw .le. 0 ) then
	nw = nw + 1
	name_w( nw ) = name
	iw = nw
      end if
* scope
	call word1( file, l1, l2 )
	name = file(l1:l2)
	file = file(l2+1:)
	is = 0
      do j = 1, ns
	if( name .eq. name_s( j ) ) is = j
	if( is .ne. 0 ) exit
      end do
      if( is .le. 0 ) then
	ns = ns + 1
	name_s( ns ) = name
* scope colour index
	k = ns
	if( k .ge. kyellow ) k = k + 1 
	kolor_s( ns ) = k
	is = ns
      end if
	write(*,*) 'Filter', iw, nw, name_w( iw )
	write(*,*) ' Scope', is, ns, name_s( is )

* lightcurve data
	call word1( file, l1, l2 )
	file = file(l1:l2)
	nt = maxp
	call loadlc3( file, 3, nt, hjd_t, dat_t, sig_t )
	write(*,*) 'File ', file(:len1(file))
	write(*,*) 'Times ', nt
	jdiv = 10000
	jd0 = 2450000
      do it = 1, nt
	hjd = hjd_t( it )
	jd = jdiv * int( hjd / jdiv )
	t = hjd - jd
	to_t( it ) = t
      end do
* report
	write(*,*) '  HJD ', hjd_t(1), hjd_t(nt)
	write(*,*) '    t ', to_t(1), to_t(nt)
	call mnmx( nt, dat_t, dmn, dmx )
	call mnmx( nt, sig_t, smn, smx )
	write(*,*) ' data', dmn, dmx	
	write(*,*) 'sigma', smn, smx	

      do it = 1, nt
* reap
	hjd = hjd_t( it )
	   t = to_t( it )
	dat = dat_t( it )
	sig = sig_t( it )
* stow
      if( no .lt. maxo ) then
	no = no + 1
	 iw_o( no ) = iw
	 is_o( no ) = is
	hjd_o( no ) = hjd
	  t_o( no ) = t
	dat_o( no ) = dat
	sig_o( no ) = sig
      else
	write(*,*) '** OBSERVATION ARRAYS FULL', no, ' :('
	cycle
      end if
      end do
	write(*,*) 'New data', nt, ' total', no, ' max', maxo
	infile_sw( is, iw ) = file

* next dataset
      end do
* report
	write(*,*)
	write(*,*) 'Waves', nw
	write(*,*) ( name_w(i)(:len1(name_w(i))), ' ', i=1,nw )
	write(*,*) 'Scopes', ns
	write(*,*) ( name_s(i)(:len1(name_s(i))), ' ', i=1,ns )
	write(*,*) ' Data', no
	call mnmx( no, t_o, tmn, tmx )
	write(*,*) ' Times', tmn, tmx
	write(*,*)

	goto 30



	goto 10

* ------------------------------------------------
* Calibration Parameter Summary
  300	write(*,*) '----------------------------------------'
	write(*,*) 'Observations', no
	call mnmx( no, t_o, tmn, tmx )
	write(*,*) 'Times', tmn, tmx
	write(*,*) 'Filters', nw,
     &		( ' ', name_w(i)(:len1(name_w(i))), i = 1, nw )
	write(*,*) ' Scopes', ns,
     &		( ' ', name_s(i)(:len1(name_s(i))), i = 1, ns )
	write(*,*) 'ROA Shape ', shape(:len1(shape)), ' Half-Width', hw
     &	, ' dt', dt
* wavelengths
      do iw = 1, nw
	write(*,*) 'Filter', iw, name_w( iw ) (:len1(name_w( iw )))
* calibration parameters
	kscl = 1
	kadd = 2
	krms = 3
      do k = 1, 3
	write(*,*) '----------------------------------------'
* scopes
      do is = 1, ns
      if( k .eq. kscl ) then
	name = 'scl'
	val = scl_sw( is, iw )
	sig = sigscl_sw( is, iw )
      else if( k .eq. kadd ) then
	name = 'add'
	val = add_sw( is, iw )
	sig = sigadd_sw( is, iw )
      else if( k .eq. krms ) then
	name = 'rms'
	val = rms_sw( is, iw )
	sig = sigrms_sw( is, iw )
      end if
	h = 'Filter '
	call append( h, name_w( iw ) )
	call append( h, 'Scope' )
	call append( h, name_s( is ) )
	call append( h, name )
	call append_datsig( h, val, sig, 3 )
	write(*,*) h(:len1(h))
      IF( .FALSE. ) THEN
	write(*,*) 'Filter ', name_w( iw ) (:len1(name_w(iw))),
     &		' Scope ', name_s( is )(:len1(name_s( is ))),
     &	' ', name(:len1(name)), val, ' +/-', sig
      END IF
* next scope is
      end do
* next parameter k
      end do
* next wavelength iw
      end do

	goto 11

* ------------------------------------------------
* Write to output files 
  400	write(*,*)

	call quicksort( no, t_o, i_o )
      do iw = 1, nw
* gather
	n = 0
      do is = 1, ns
	nt_s( is ) = 0
      end do
      do i = 1, no
	io = i_o( i )
      if( iw_o( io ) .eq. iw ) then
	hjd = hjd_o( io )
	  t =   t_o( io )
	dat = dat_o( io )
	sig = sig_o( io )

	cal = cal_o( io )
	sigcal = sigcal_o( io )
	rmscal = rmscal_o( io )
	sigoldcal = sigoldcal_o( io )
	is = is_o( io )
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
* stow
	n = n + 1
	hjd_t( n ) = hjd
	  t_t( n ) = t
	dat_t( n ) = cal
	sig_t( n ) = sigcal
	plot1( n ) = rms
	plot2( n ) = is
	plot3( n ) = add
	plot4( n ) = scl
	plot5( n ) = sigoldcal
	nt_s( is ) = nt_s( is ) + 1
      end if
      end do


      if( n .gt. 0 ) then

* prefix
	l = len1( obsfile )
	write(*,*) 'Input File: ', obsfile(:l)
      if( l .gt. 0 ) then
	prefix = obsfile
	l = index( prefix, '.' )
	if( l .ge. 2 ) prefix = prefix(:l-1)
      end if
	call inq( 'Output file PREFIX:', prefix, prefix)
	call word1( reply, l1, l2 )
	reply = reply( l1 : l2 )

* output file
	file = ' '
	l = len1( prefix )
	if( l .gt. 0 ) file = prefix(:l) // '_'
	name = name_w( iw )
	file = file(:len1(file)) // name(:len1(name)) // '_dream.dat'

	call inq( 'Output file:', file, file )
	call word1( file, l1, l2 )
	file = file( l1 : l2 )
	reply = file
	call upcase( reply, reply )

      if( file .ne. ' ' .and. file .ne. 'NONE' ) then
	
* header
	h_h( 1 ) = 'File: ' // file(:len1(file))
	h_h( 2 ) = '  Col 1 : HJD'
	h_h( 3 ) = '  Col 2 : flux'
	h_h( 4 ) = '  Col 3 : sqrt( col4^2 + col5^2 )'
	h_h( 5 ) = '  Col 4 : rms statistical error'
	h_h( 6 ) = '  Col 5 : rms systematic error'
	h_h( 7 ) = '  Col 6 : scope (see key below)'
	h_h( 8 ) = '  Col 7 : add'
	h_h( 9 ) = '  Col 8 : scl'
	nh = 9
* version
	h = 'DREAM'
	call append( h, 'version' )
	call append( h, version )
	call append( h, 'iteration' )
	call append_i( h, iter )
	nh = nh + 1
	h_h( nh ) = h
* lightcurve model
	h = 'ROA shape ' // shape(:len1(shape))
	call append( h, 'half-width' )
	call append_n( h, hw, 5 )
	call append( h, 'dt' )
	call append_n( h, dt, 5 )
	nh = nh + 1
	h_h( nh ) = h

* band
	call avgrmsx( n, dat_t, sig_t, good, avg, rms, sigavg, sigrms )
	h = 'Band '
	call append( h, name_w( iw ) )
	call append( h, 'Nt' )
	call append_i( h, n )
	call append( h, 'avg' )
	call append_datsig( h, avg, sigavg, 3 )
	call append( h, 'rms' )
	call append_datsig( h, rms, sigrms, 3 )
	nh = nh + 1
	h_h( nh ) = h
      do is = 1, ns
	add = add_sw( is, iw )
	scl = scl_sw( is, iw )
	rms = rms_sw( is, iw )
	sigadd = sigadd_sw( is, iw )
	sigscl = sigscl_sw( is, iw )
	sigrms = sigrms_sw( is, iw )

	h = 'Scope'
	call append_i( h, is )
	call append( h, name_s( is ) )

	call append( h, 'Nt' )
	call append_i( h, nt_s( is ) )

	call append( h, 'add' )
	call append_datsig( h, add, sigadd, 3 )

	call append( h, 'scl' )
	call append_datsig( h, scl, sigscl, 3 )

	call append( h, 'rms' )
	call append_datsig( h, rms, sigrms, 3 )

	nh = nh + 1
	h_h( nh ) = h
      end do
	nh = min( nh, maxh )

* report
	write(*,*)
	write(*,*) 'Headlines:'
      do i = 1, nh
	h = h_h( i )
	write(*,*) ' ', h( : len1(h) )
      end do
	write(*,*)
	write(*,*) 'First data rows:'
      do i = 1, min( n, 5 )
	write(*,*) hjd_t(i), dat_t(i), sig_t(i)
     &	, plot5(i), plot1(i), plot2(i), plot3(i), plot4(i)
      end do
	write(*,*)

* write to disk
	write(*,*) 'dumplc8h(', nh, 8, n, ' )'
	call dumplc8h( file, nh, h_h, 8, n, hjd_t, dat_t, sig_t
     &	, plot5, plot1, plot2, plot3, plot4 )
      end if
      end if
* next wavelength iw
      end do

	goto 11


	end
**************************************************
	subroutine fitmin( n, x, y, i0, x0, y0, b, c )
* find minimum and fit parabola to 3 lowest points
*   y = y0 + b (x-x0) + c (x-x0)^2
* In:	n	i4 number of data
*	x(n)	r4 x values (increasing)
*	y(n)	r4 y values 
* Out:	i0	i4 index of lowest y
*	x0	r4 x value at min (inside grid)
*	y0	r4 y value at min
*	b	r4 slope
*	c	r4 curvature (0 if lowest y at grid edge)
* 2020 Aug Keith Horne @ Crail
c	implicit real*8 (a-h, o-z)
	dimension x( n ), y( n )
* no data
	x0 = 0.d0
	y0 = 0.d0
	c = 0.d0
	if( n .le. 0 ) return
* find minimum
	x0 = x( 1 )
	i0 = 1
      do i = 1, n
	if( y( i ) .lt. y( i0 ) ) i0 = i
      end do 	
	x0 = x( i0 )
	y0 = y( i0 )
	write(*,*) 'Min at grid', i0, ' of', n
	write(*,*) 'Min on grid x', x0, ' y', y0
	c = 0.d0
* no data for curvature
	if( n .lt. 3 ) return
* domain boundary
	ip = i0 + 1
      if( i0 .eq. 1 ) then
	b = ( y( ip ) - y0 ) / ( x( ip ) - x0 )
	return
      end if
	im = i0 - 1
      if( i0 .eq. n ) then
	b = ( y0 - y( im ) ) / ( x0 - ( im ) )
	return
      end if
* quadratic fit
	xm = x( im ) - x0
	xp = x( ip ) - x0
	ym = y( im ) - y0
	yp = y( ip ) - y0
	write(*,*) im, xm, ym
	write(*,*) ip, xp, yp
	xm2 = xm * xm
	xp2 = xp * xp
	det = xm * xp2 - xp * xm2
	b = ( xp2 * ym - xm2 * yp ) / det
	c = ( xm * yp - xp * ym ) / det
	dx = - b / ( c + c )
	x0 = x0 + dx
	y0 = y0 + ( b + c * dx ) * dx
	write(*,*) 'Min of fit x', x0, ' y', y0
	write(*,*) 'Slope', b, ' Curvature', c
	return
	end
