* simulate dwarf nova lightcurve for Kronos

* 2001 Sep Keith Horne @ St.Andrews
	parameter ( maxt = 100000 )

	real*4 t( maxt )
	real*4 phase( maxt )
	real*4 cycle( maxt )
	real*4 flux( maxt )

	real*4 flicker( maxt )

	real*4 pic( maxt )
	real*4 pic1( maxt )
	real*4 pic2( maxt )
	real*4 pic3( maxt )
	real*4 pic4( maxt )

	real*8 hjd( maxt ), time0
	real*8 period, hjd0, e, dt
	character*75 file
	character*75 xlabel, ylabel, title
	logical pgagain

* constants
	hour = 3600.
	day = hour * 24.

* binary parameters
	phr = 3.9146
	period = phr * hour / day
* z cha parameters
	pmin = 107.
* fluxes
	fwd = 0.7
	fbs = 1.1
	fq = 0.7
	fo = 28.
	fk = 0.2

	period = pmin * 60. / day
	time0 = 0.
	hjd0 = time0 + 0.5 * period

* superhump
	psh = 1.043784 * period
	fsh = 0.3
	write(*,*) 'period:', phr, period

	totaldays = 20.
	dt = ( hour / 60. ) / day
	nbin = 100
	dt = period / nbin
	nt = 1 + ( totaldays / dt )
	nt = min( nt, maxt )
	nt = nbin * ( nt / nbin )
	write(*,*) 'nt:', nt, maxt


* flickering
	iseed = 47
	tfast = 5. / day / dt
	tslow = 0.2 * period / dt
	frms = 0.3
	write(*,*) 'flicker timescales:', tfast, tslow, ' rms:', frms
      do i=1,nt
	flicker(i) = ranw( iseed, tfast )
	pic(i) = flicker(i)
      end do
	call mnmx( nt, flicker, fmn, fmx )
	write(*,*) 'fmn,fmx:', fmn, fmx
	nbox = max( 1, tslow / dt )
	call boxfilt( nt, pic, nbox )
      do i=1,nt
	flicker(i) = flicker(i) - pic(i)
      end do
	call mnmx( nt, flicker, fmn, fmx )
	write(*,*) 'fmn,fmx:', fmn, fmx
	call avgrms( nt, flicker, avg, rms )
	write(*,*) 'flicker:', avg, rms
      do i=1,nt
	eta = ( flicker(i) - avg ) / rms
	flicker(i) = exp( frms * eta )
      end do
	call mnmx( nt, flicker, fmn, fmx )
	write(*,*) 'fmn,fmx:', fmn, fmx

      do i=1,nt
	hjd(i) = time0 + i * dt
	t(i) = hjd(i)

* superhump
	e = ( hjd(i) - hjd0 ) / psh
	phish = e - nint( e )
* binary
	e = ( hjd(i) - hjd0 ) / period
	phi = e - nint( e )

	oblc = dnob( t(i) ) * exp( fsh * shlc( phish ) )
	fdisk = ( fq + (fo-fq) * oblc ) * disklc( phi ) 
	fdisk = fdisk * exp( 0.2 * alog( flicker(i) ) )

	flx = fk + fdisk + fwd * wdlc( phi ) 
     &	+ fbs * flicker(i) * bslc( phi )

	flux(i) = flx
	phase(i) = e
* next time
      end do


* normalize cycles

	ncycle = nt / nbin
      do i=1,ncycle
	cycle(i) = i
	j = 1 + (i-1) * nbin
	call peqp( nbin, pic(j), flux(j) )
* avg and rms
	call avgrms( nbin, pic(j), avg, rms )
	write(*,*) i, avg, rms
* divide by mean
	call peqpmc( nbin, pic(j), pic(j), 1./avg )
      end do
	call mnmx( nt, pic, pmn, pmx )
	write(*,*) 'normalized lightcurves: ', ncycle, pmn, pmx

* kronos sampling
      do i=1,nt
	pic1(i) = pic(i) * viskronos( t(i) )
      end do
* hst sampling
      do i=1,nt
	pic2(i) = pic(i) * vishst( t(i) )
      end do
* ground sampling
      do i=1,nt
	pic3(i) = pic(i) * visground( t(i) )
      end do
* chandra sampling
      do i=1,nt
	pic4(i) = pic(i) * vischandra( t(i) )
      end do

* plot lightcurves ------------------------------

	nloop = 4
      if( nloop .eq. 1 ) then
	nt = maxt
	file = 'cvlc.dat'
	call load1d2( file, 2, nt, t, ab )
	write(*,*) 'times:', nt, t(1), t(nt)
      end if

  100	call pgstart( 1, nloop, new )
	csize = 1.5
	lwide = 2
	call pgsch( csize )
	call pgslw( lwide )
      do loop=1,nloop
      if( loop .eq. 1 ) then
	i = 1
	np = nt
      else if( loop .eq. 2 ) then
	want = t(1) + period
	call lookup( nt, t, want, i1, i2, part )
	i = 1
	np = 10.5 * i2
      else if( loop .eq. 3 ) then
	call lookup( nt, t, 2.3, i1, i, part )
      else if( loop .eq. 4 ) then
	call lookup( nt, t, 7., i1, i, part )
      end if
	call mnmx( np, t(i), xmn, xmx )
	call mnmx( np, flux(i), ymn, ymx )
	write(*,*) xmn, xmx, ymn, ymx
	ymn = 0.
	call pgsci( 5 )
	call pgenv( xmn, xmx, ymn, ymx, 0, 0 )
	call pgsci( 7 )
	xlabel = 'time (d)'
	ylabel = 'f\d\gn\u (mJy)'
	title = 'Simulated Dwarf Nova Lightcurve (Z Cha)'
	call pglabel( xlabel, ylabel, title )
	call pgsci( 1 )
	call pgline( np, t(i), flux(i) )
      end do
	call pgstop
	if( pgagain() ) goto 100


* greyscale plot ----------------------------------------

* format
	kronos = 1
	ihst = 2
	ichandra = 3
	iground = 4
	nloop = 4

	nx = nloop
	ny = 1

* open
  200	call pgstart( nx, ny, new )
	csize = 1.5
	lwide = 2
	call pgsch( csize )
	call pgslw( lwide )

      do loop=1,nloop

      if( loop .eq. kronos ) then
	call peqp( nt, pic, pic1 )
	title = 'KRONOS'
      else if( loop .eq. ihst ) then
	call peqp( nt, pic, pic2 )
	title = 'HST'
      else if( loop .eq. iground ) then
	call peqp( nt, pic, pic3 )
	title = 'ground'
      else if( loop .eq. ichandra ) then
	call peqp( nt, pic, pic4 )
	title = 'Chandra'
      end if

* window
	ymn = cycle( 1 )
	ymx = cycle( ncycle )
	xmn = phase(1) - 0.5 / nbin
	xmx = xmn + 1.
	write(*,*) xmn, xmx, ymn, ymx
	call pgsci( 5 )
	call pgenv( xmn, xmx, ymn, ymx, 0, 0 )
* label
	xlabel = 'binary phase'
	ylabel = 'binary cycle'
	
	call pgsci( 7 )
	call pglabel( xlabel, ylabel, title )
* greyscale
	bg = 0.
	vg = pmx * 1.2
	write(*,*) 'fg,bg:', fg, bg
	write(*,*) nbin, phase(1), phase(nbin)
	write(*,*) ncycle, cycle(1), cycle(ncycle)

	call pgpic( 'G', pic, nbin, ncycle, 1,nbin, 1,ncycle,
     &		fg, bg, phase, cycle )
	call pgsci( 5 )
	call pgbox( 'bcst', 0., 0, 'bcst', 0., 0 )

* next loop
      end do

	call pgstop
	if( pgagain() ) goto 200


	write(*,'(A,$)') ' Output file : '
	read(*,'(A)') file
	if( file.ne. ' ' ) call dump1d2( file, 2, nt, t, flux )
	stop
	end

**********************************
	function dnob( t )
* white dwarf lightcurve (0-1)
	logical firstcall /.true./
	parameter ( maxt = 100 )
	real*4 tlc( maxt ), flc( maxt )
      if( firstcall ) then
	tlc(1) = -9999.
	flc(1) = 0.
	tlc(2) = 2.374
	flc(2) = 0.
	tlc(3) = 3.74
	flc(3) = 0.8
	tlc(4) = 5.709
	flc(4) = 1.
	tlc(5) = 11.749
	flc(5) = 0.8
	tlc(6) = 13.749
	flc(6) = 0.7
	tlc(7) = 17.317
	flc(7) = 0.
	tlc(8) = 99999.
	flc(8) = 0.
	nt = 8
	firstcall = .false.
      end if

	call lookup( nt, tlc, t, i1, i2, part )
	flx = flc(i1) * (1.-part) + flc(i2) * part
	dnob = flx

	return
	end
**********************************
	function vishst( t )
* hst visibility
	logical firstcall /.true./
      if( firstcall ) then
	day = 3600. * 24.
	pminutes = 96.
	phst = pminutes * 60. / day
	duty = 59. / pminutes
	t0 = ran( iseed ) * phst
	firstcall = .false.
      end if
	phi = ( t - t0 ) / phst
	phi = phi - nint( phi )
	vis = 1.
	if( abs( phi ) .gt. duty/2. ) vis = 0.
	vishst = vis
	return
	end
**********************************
	function vischandra( t )
* chandra visibility
	logical firstcall /.true./
      if( firstcall ) then
	phr = 47.3 + 16.7
	pday = phr / 24.
	duty = 47.35 / phr
	t0 = ran( iseed ) * pday
	firstcall = .false.
      end if
	phi = ( t - t0 ) / pday
	phi = phi - nint( phi )
	vis = 1.
	if( abs( phi ) .gt. duty/2. ) vis = 0.
	vischandra = vis
	return
	end
**********************************
	function visground( t )
* ground visibility
	logical firstcall /.true./
      if( firstcall ) then
	duty = 10./24.
	firstcall = .false.
      end if
	phi = t - nint(t)
	vis = 1.
	if( abs( phi ) .gt. duty/2. ) vis = 0.
	visground = vis
	return
	end
**********************************
	function viskronos( t )
* kronos visibility
	logical firstcall /.true./
      if( firstcall ) then
	duty = 1.
	pday = 13.7
	duty = 1. - ( 3. / 24. / pday )
	t0 = pday * ran( iseed )
	firstcall = .false.
      end if
	phi = ( t - t0 ) / pday
	phi = phi - nint(phi)
	vis = 1.
	if( abs( phi ) .gt. duty/2. ) vis = 0.
	viskronos = vis
	return
	end
**********************************
	function wdlc( phi )
* white dwarf lightcurve (0-1)
	logical firstcall /.true./
      if( firstcall ) then
	call wdcontacts( w1, w2, w3, w4 )
	dwe = w4 - w3
	firstcall = .false.
      end if
	vis = ( abs( phi ) - w3 ) / dwe
	wdlc = max( 0., min( 1., vis ) )
	return
	end
**********************************
	subroutine wdcontacts( w1, w2, w3, w4 )
* z cha quiescence from Wood et al
	dphi = 0.0535
	dw = 0.0085
	w1 = - ( dphi +  dw ) /2.
	w2 = w1 + dw
	w3 = - w2
	w4 = - w1
	return
	end
**********************************
	subroutine bscontacts( b1, b2, b3, b4 )
* z cha quiescence from Wood et al
	dphi = 0.0902
	bi = -0.012
	be =  0.078
	dbi = 0.0104
	dbe = 0.0089
	b1 = bi - dbi/2.
	b2 = b1 + dbi
	b3 = be - dbe/2.
	b4 = b3 + dbe
	return
	end
**********************************
	function disklc( phi )
* disk dwarf lightcurve (0-1)
	logical firstcall /.true./
      if( firstcall ) then
	call wdcontacts( w1, w2, w3, w4 )
	sig = ( w3 + w4 ) / 2.
	deep = 0.8
	firstcall = .false.
      end if
	xx = phi / sig
	g = exp( -0.5 * xx**2 )
	p = 0.05
	vis = max( 0., g - p ) / ( 1. - p )
	disklc = 1. - deep * vis
	return
	end
**********************************
	function bslc( phi )
* brightspot lightcurve normalized to 1 at peak
	logical firstcall /.true./
      if( firstcall ) then
	twopi = 8. * atan2( 1., 1. )
	call bscontacts( b1, b2, b3, b4 )
	dbi = b2 - b1
	dbe = b4 - b3
	phihump = -0.13
	fiso = 0.3
	fhump = 1. - fiso
	firstcall = .false.
      end if
	vis1 = ( b2 - phi ) / dbi
	vis2 = ( phi - b3 ) / dbe
	vis = max( vis1, vis2 )
	bslc = 0.
	if( vis .le. 0. ) return
	vis = min( 1., vis )
	radian = twopi * ( phi - phihump )
	flx = fiso + fhump * max( 0., cos( radian ) ) ** 1.3
	bslc = flx
	return
	end
**********************
	function shlc( phi )
	twopi = 8. * atan2( 1., 1. )
	radian = twopi * phi
	c1 = ( 1 + cos( radian ) ) / 2.
	c2 = ( 1 + cos( radian - 9. ) ) / 2.
	c3 = ( 1 + cos( radian + 3. ) ) / 2.
	shlc = c1 ** 3 - 0.1 * c2**7 + 0.2 * c3**8
	return
	end
