* plens -- micro-lens lightcurve analysis program
* 1998 Aug Keith Horne @ St-Andrews
* 2002 Nov KDH @ StAnd - blending, param editor
* 2003 Jul KDH @ StAnd - multi-site
* 2003 Sep KDH @ StAnd - prior(f)
* 2004 Jun KDH @ StAnd - prior(te,A0)
* 2006 Jan KDH @ StAnd - scalerb_l(l)
* 2007 Mar KDH @ Goettingen - guesslens more robust
* 2007 Mar KDH @ Goettingen - syserr
* 2007 Mar KDH @ StAnd - chiscale->erbscl
* 2007 Mar KDH @ StAnd - parameter indices kt0,kte,...
* 2007 Oct KDH @ Christchurch -- Chi2(t,f)
* 2008 Feb KDH @ StAnd - call die rather than pause
* 2008 Jul KDH @ StAnd - fit fluxes rather than mags
* 2008 Jul KDH @ SAAO - Priorities on mag plot
* 2008 Jul KDH @ SAAO - blend and sigma grids working.
* 2008 Jul KDH @ SAAO - parameter covariances
* 2008 Jul KDH @ SAAO - parameter output file
* 2008 Aug KDH @ StA - menus in gridsigma and gridblend
* 2009 May KDH @ StA - spawn command
* 2009 May KDH @ StA - baseline rebinning
* 2009 May KDH @ StA - split-gaussian prior on b=f_b/f_s
* 2009 Jun KDH @ StA - white background and site-specific colours
* 2009 Jul KDH @ StA - time zone
* 2009 Jul KDH @ StA - sort hjd on input if required
* 2009 Jul KDH @ StA - multi-site blend grid, 
* 2009 Nov KDH @ StA - move subroutines to plens.a
* 2009 Dec KDH @ StA - gridblend - dirty bug fix, diagnostic plot
* 2010 Jun KDH @ StA - use gmtime for UT
* 2010 Jun KDH @ StA - date_and_time replaces idate (unknown to ifort)

	include 'plens.inc'
* version
	version = '2009-Dec-01'
	version = '2010-Mar-19'
	version = '2010-Mar-24'
	version = '2010-Apr-23'
	version = '2010-May-07'
	version = '2010-Jun-08'
	version = '2011-Mar-23'
	version = '2011-Jun-20'
	version = '2011-Jun-22'
* mute screen messages
	muteamoeba = 1
	muteamoebaplot = 1
	mutefit0 = 1

* intro
	write(*,*) '**************************************************'
	write(*,*) '* PLENS - fit binary lens model to lightcurves   *'
	write(*,*) '*       - assess lens star planet detect zones   *'
	write(*,*) '* 1998 Aug - Keith Horne @ St.Andrews            *'
	write(*,*) '* 2002 Nov - blending                            *'
	write(*,*) '* 2003 Jul - multi-site baseline mags            *'
	write(*,*) '* 2003 Sep - Prior(f)                            *'
	write(*,*) '* 2003 Dec - seeing                              *'
	write(*,*) '* 2004 Jun - Prior(tE,A0)                        *'
	write(*,*) '* 2007 Mar - more robust guesslens               *'
	write(*,*) '* 2007 Oct - Chi2(x,y) (lnu,theta) (t,f) maps    *'
	write(*,*) '* 2008 Jul - fit fluxes, Priority on mag plot    *'
	write(*,*) '* 2008 Jul - parameter covariancess and writeout *'
	write(*,*) '* 2009 May - baseline rebin and blend prior      *'
	write(*,*) '* ', version(:len1(version)), ' = current version' 
	write(*,*) '**************************************************'

* time zone
	timezone = 1.

* initialize
	call initplens

* pre-calculate the 4,5,6-param chi^2 values
   10 if( ndata .gt. 0 ) then
	title = ' '
	call title_chi2( title )
	call nogreek( title )
	call noupdn( title )
	call badness( bof, nd, np )
      end if
* menu
   	write(*,*) 'PLENS menu: ====================================='
	call system( 'date' )
	call now( iyr, imon, iday, ihr, imin, isec )
	write(*,'(a,i4,5(a,i2.2))')
     &		'now ', iyr, '-', imon, '-', iday
     &		,' UT ', ihr, ':', imin, ':', isec

* data
	title = ' D ... Data:'
	call append_i( title, ndata )
      if( ndata .gt. 0 ) then
	call append( title, datfile(:len1(datfile)) )
      if( event .ne. ' ' ) then
	call append( title, event(:len1(event)) )
      end if
      end if
	write(*,*) title(:len1(title))

* multi-site report
	call multisitereport

* params
	write(*,*) ' E ... Edit parameters'
	call showpar
* time zone
      if( .false. ) then
	call now( iy, im, id, ith, itm, its )
	tlocal = ih + ( im + is/60.d0)/60.d0
c	call get_dmy( id, im, iy )
c	call get_hmst( ith, itm, its, tlocal )
	utnow = tlocal - timezone
	write(*,*) ' TZ .. Time Zone:', timezone
     &		, ' local', tlocal, ' UT', utnow
      end if
* time range
	t0 = hjd0 - mjd0
	write(*,*) ' T ... Time range:', tbeg, tend
     &	, ' (t-t0)/te', (tbeg-t0)/te, (tend-t0)/te

* fits
      if( ndata .gt. 0 ) then
	write(*,*) ' X ... Bin baseline', binday, ' days.'
	title = ' '
	call title_chi2( title )
	call nogreek( title )
	call noupdn( title )
	call noeqbl( title )
	write(*,*) ' F ... Fit amoeba ', title(:len1(title))
	write(*,*) ' W ... Write parameter file'
	write(*,*) ' WXY,WUA,WTF write Chi^2 map files'
	title = ' '
	call title_sig0( title )
	call title_f( title )
	call title_pr_f( title )
	call noeqbl( title )
	write(*,*) ' S ... Sigma grid ', title(:len1(title))
	title = ' '
	call title_b( title )
	call noeqbl( title )
	write(*,*) ' B ... Blending (BB=grid) ', title(:len1(title))
      end if
* plots
* magnification
	title = ' '
	call title_anow( title )
	call title_a0pm( title )
	write(*,*) ' A ... A(t) plot ', title(:len1(title))
* magnitude
	title = ' '
	call title_priority( title )
	call noupdn( title )
	write(*,*) ' M ... Mag(t) ', title(:len1(title))

      if( ndata .gt. 0 ) then
	write(*,*) ' R ... Residuals'   
	write(*,*) ' N ... Normalized residuals'
	write(*,*) ' G ... Gaussian residuals test'
	write(*,*) ' Y ... Sky correlations'
	write(*,*) ' Z ... Seeing correlations'
* chi2(x,y)
c	write(*,*) ' C ... Chi2(x,y|q) ', nxp, nyp, qxy, chi2xy, nint(cpuxy)
	title = ' C ... Chi2(x,y|q)'
	call append_i( title, nxp )
	call append_i( title, nyp )
	call append( title, 'q' )
	call append_n( title, qxy, 3 )
	call append( title, 'Chi2' )
	call append_f( title, chi2xy, '(f10.2)' )
	call append( title, 'cpu' )
	call append_i( title, nint(cpuxy) )
	write(*,*) title(:len1(title))
* chi2(u,a)
c	write(*,*) ' L ... Chi2(u,a|q)', nup, nap, qua, chi2ua, nint(cpuua)
	title = ' L ... Chi2(u,a|q)'
	call append_i( title, nup )
	call append_i( title, nap )
	call append( title, 'q' )
	call append_n( title, qua, 3 )
	call append( title, 'Chi2' )
	call append_f( title, chi2ua, '(f10.2)' )
	call append( title, 'cpu' )
	call append_i( title, nint(cpuua) )
	write(*,*) title(:len1(title))
* chi2(t,f)
c	write(*,*) ' V ... Chi2(t,f|q) ', ntp, nfp, qtf, chi2tf, nint(cputf)
	title = ' V ... Chi2(t,f|q)'
	call append_i( title, ntp )
	call append_i( title, nfp )
	call append( title, 'q' )
	call append_n( title, qtf, 3 )
	call append( title, 'Chi2' )
	call append_f( title, chi2tf, '(f10.2)' )
	call append( title, 'cpu' )
	call append_i( title, nint(cputf) )
	write(*,*) title(:len1(title))
* chi2(f,q)
c	write(*,*) ' K ... Chi2(f,q)', nf, nqp, qfq, chi2fq, nint(cpufq)
	title = ' K ... Chi2(f,q)'
	call append_i( title, nf )
	call append_i( title, nqp )
	call append( title, 'q' )
	call append_n( title, qfq, 3 )
	call append( title, 'Chi2' )
	call append_f( title, chi2fq, '(f10.2)' )
	call append( title, 'cpu' )
	call append_i( title, nint(cpufq) )
	write(*,*) title(:len1(title))

	write(*,*) ' P ... P(det|a)'
	write(*,*) ' H ... Hone binary-lens fits'
	write(*,*) ' U ... U-q grid '
      end if
	write(*,*) ' I ... Image positions'
	write(*,*) ' O ... Orbit P(x,y|a)'
* fini
	write(*,*) ' $ ... Spawn a command to the operating system'
	write(*,*) ' Q ... Quit PLENS'

* command node
   15	call inq( 'PLENS> ', ' ', reply )

* spawn a command to operating system
	call word1( reply, l1, l2 )
      if( l1 .ge. 1 ) then
      if( reply(l1:l1) .eq. '$' ) then
	reply = reply(l1+1:)
	call spawn( reply(:len1(reply)) )
	goto 15
      end if
      end if

	call upcase( reply, reply )
	if( reply .eq. ' ' ) goto 10

* change mute options

	call word1( reply, l1, l2 )
      if( reply(l1:l1+3) .eq. 'MUTE' ) then
   20	write(*,*) 'Mute diagnostics:'
	write(*,*) 'F ... mute Fit0', mutefit0
	write(*,*) 'A ... mute Amoeba', muteamoeba
	write(*,*) 'P ... mute amoebaPlot', muteamoebaplot
	reply = ' '
	call inq( 'Which to change', reply, reply )
	call upcase( reply, reply )
      if( reply .eq. 'F' ) then
	mutefit0 = max( 0, 1 - mutefit0 )
      else if( reply .eq. 'A' ) then
	muteamoeba = max( 0, 1 - muteamoeba )	
      else if( reply .eq. 'P' ) then
	muteamoebaplot = max( 0, 1 - muteamoebaplot )	
      else
	goto 15
      end if
	goto 20
      end if

* quit
      if( reply(1:1) .eq. 'Q' ) then
	reply = 'N'
	call inq( 'Do you REALLY want to quit PLENS ?', reply, reply )
	call upcase( reply, reply )
	if( reply .eq. 'Y' ) stop 'PLENS retired.'

* timezone
      else if( reply(1:2) .eq. 'TZ' ) then

* local date/time (computer clock)
	write(*,*) 'local computer date/time:'
	call system( 'date' )
	call now( iyr, imon, iday, ihr, imin, isec )
	write(*,'(a,i4,5(a,i2.2))')
     &		'now ', iyr, '-', imon, '-', iday
     &		,' UT ', ihr, ':', imin, ':', isec

c	tlocal = ih + ( im + (is + ims/1000.d0)/60.d0)/60.d0
c	utnow = tlocal - timezone
c	write(*,*) 'local', tlocal, ' UT', utnow, ' timezone', timezone
	write(*,*) 'BST +1  GMT 0  EST -5  CST -6  MST -7  PST -8'
	ival = 0
	call inqi( 'Time Zone (hours ahead of UTC):', ival, ival )
	timezone = ival

* plots
      else if( reply(1:1).eq. 'A' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'I' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'M' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'R' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'N' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'P' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'G' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'Y' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'Z' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'V' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'C' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'K' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'L' ) then
	call mplot( reply )
      else if( reply(1:1) .eq. 'O' ) then
	call mplot( reply )
      else if( reply(1:2) .eq. 'XY' ) then
	call mplot( reply )
      else if( reply(1:2) .eq. 'UA' ) then
	call mplot( reply )
      else if( reply(1:2) .eq. 'TF' ) then
	call mplot( reply )
      else if( reply(1:2) .eq. 'FQ' ) then
	call mplot( reply )
* rebin
      else if( reply(1:2) .eq. 'X' ) then
	call rebindata
* params
      else if( reply(1:1) .eq. 'E' ) then
	call editparms
c	npfit = 1
	call needs
* time range
      else if( reply(1:1) .eq. 'T' ) then
	t0 = hjd0 - mjd0
	write(*,*) 'HJD0-MJD0 =', hjd0, ' -', mjd0, ' =', t0
	call edittimes
* data
      else if( reply(1:1) .eq. 'D' ) then
	call newdata
	call needs
* fits
      else if( ndata .gt. 0 ) then
      if( reply(1:1) .eq. 'F' ) then
	npf = 0
	call fitlens( npf )
	call needs
* covariances
	call covcalc
	call showbigcov( 0.5 )
	call showcov
	call showparsig
	goto 15
* error bars
      else if( reply(1:1) .eq. 'S' ) then
	call gridsigma
	call needs
* blending
      else if( reply(1:1) .eq. 'B' ) then
      if( reply(1:2) .eq. 'BB' ) then
	call gridblend
	call needs
      else
	call mplot( reply )
      end if

* u-q grid
      else if( reply(1:1) .eq. 'U' ) then
	call griduq
	call needs

* hone binary lens fits
      else if( reply(1:1) .eq. 'H' ) then
	call hone
	call needs

* write params to file
      else if( reply(1:1) .eq. 'W' ) then
	call writeout( reply )
	call needs
	goto 15
      end if
      else
	goto 10
      end if
	goto 10
	end
