* 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