; IDL routines used in the 30210 timing analysis exercise
;
; Søren Brandt 2014 
;
; return time: IJD (INTEGRAL Julian Day) in units of days
;        t: time in seconds
; (x,y,e) is coordinate and pulse height in JEM-X
; im is a 2D histogram of (x,y), the detector image
; of events from Crab observation in orbit 40
pro step2,name,time,t,x,y,e,im
  crab40files,2,res
  getevdata,res(1),2,time,x,y,e
  im=intarr(256,256)
  t=(time-time(0))*86400.0
  name=res(2)  
  for i=0L,n_elements(x)-1 do im(x(i),y(i))+=1
  print,'X-ray event data file read: ', name
  print,''
  print,'Returned values:'
  print,'time: IJD (INTEGRAL Julian Day) in units of days'
  print,'(x,y,e) is coordinate and pulse height (energy) in JEM-X'
  print,' '
  print,'im is a 2D histogram of (x,y), the detector image
  print,' '
  print,'data are from INTEGRAL revolution 40 in 2003'
end

; analyse the JEM-X data of revolution 40 (in Feb 2003)
; determines the Crab frequency in each science window after
; barycentric corrections
; INPUT: jemx  the JEM-X unit to be read, allowed is 1 or 2
; OUTPUT: 
;    tt    array containing the mean IJD of each observation
;    freq  array of the best Crab frequency for each observation (Hz)
;    dtt   array of duration of each observation
;    dfreq uncertainty of the frequency as the width of Gauss fit 
pro do40,jemx,tt,freq,dtt,dfreq
  ; find the available observations
  crab40files,jemx,res
  n=n_elements(res)
  freq=dblarr(n)
  dfreq=freq
  tt=dblarr(n)
  dtt=tt
  ; loop over observations
  for i=0,n-1 do begin
    ; get event time tags of the observation
    getdata,res(i),jemx,time
    
    ; get the orbit information  
    readorbit,40,epoch,xyz,r,rdist,smaxis
    ; do the sat to barycenter corrections in the Crab direction
    sat2barycenter,xyz,epoch,epochbc,time,timebc
    ; tbc is barycenter corrected in units of seconds
    tbc=(timebc-timebc(0))*86400.0

    if (n_elements(flag) gt 0) then tbc=(time-time(0))*86400.0 
    plot,histogram(tbc)
    ; do pwoer spectrum to determine the fundamental
    powerd,tbc,100,f,pd
    plot,f,pd
    w=where(f gt 29 and f lt 30)
    print,max(pd(w),n)
    pmax=1/f(w(n))
    ; do period folding near the peak
    psearch,pmax,0.0000001,50,tbc,p,ff,h,pmax
    ; Gauss fit the merit function as function of frequency
    r=gaussfit(1.0/p,ff,nterms=3,a)
    print,a
    ; do another finer 10 ns period folding
    psearch,pmax,0.00000001,50,tbc,p,ff,h,pmax
    ; do another Gauss fit in frequncy
    r=gaussfit(1.0/p,ff,nterms=4,a,sigma=da)
    print,a
    ; save the fitted frequency
    freq(i)=a(1)
    ; save the Gauss width as the uncertainty
    dfreq(i)=da(1)
    ;dfreq(i)=a(2)
    ; save the mean time
    tt(i)=mean(time)
    ; save the duration of the observation
    dtt(i)=max(time)-min(time)
  endfor
end

; analyse rev 239
pro do239,time,tbc,h,pmax,meantbc,dmeantbc
  crab239files,res
  getdata,res(0),1,time1
  getdata,res(2),1,time
  time1=[time1,time]
  getdata,res(1),2,time2
  getdata,res(3),2,time
  time2=[time2,time]
  time=[time1,time2]
  time=time(sort(time))
  readorbit,239,epoch,xyz,r,rdist,smaxis
  sat2barycenter,xyz,epoch,epochbc,time,timebc
  meantbc=mean(timebc)
  dmeantbc=max(timebc)-min(timebc)
  tbc=(timebc-timebc(0))*86400.0
  plot,histogram(tbc)
  powerd,tbc,200,f,pd
  plot,f,pd
  w=where(f gt 29 and f lt 30)
  print,max(pd(w),n)
  pmax=1/f(w(n))
  psearch,pmax,0.000001,50,tbc,p,ff,h,pmax
  r=gaussfit(1.0/p,ff,nterms=3,a)
  print,a
  plot,p,ff,psym=10
  oplot,p,r,co=255
  psearch,pmax,0.0000001,50,tbc,p,ff,h,pmax
  r=gaussfit(1.0/p,ff,nterms=3,a)
  print,a
  plot,p,ff,psym=10
  oplot,p,r,co=255
  psearch,pmax,0.00000001,50,tbc,p,ff,h,pmax
  r=gaussfit(1.0/p,ff,nterms=4,a)
  print,a
  plot,p,ff,psym=10
  oplot,p,r,co=255
end

; return the file names of data files for rev 40
; INPUT  jemx   JEM-X unit 1 or 2
; OUTPUT res   array of file names
pro crab40files,jemx,res
  if (jemx eq 1) then begin 
    name='/home/sb/jemx/Crab/40_r3/00400002???0.001/jmx1_events.fits'
    res=findfile(name,count=n)
    name='/home/sb/jemx/Crab/40_r3/00400003???0.001/jmx1_events.fits'
    res1=findfile(name,count=n)
    res=[res,res1]
    name='/home/sb/jemx/Crab/40_r3/00400004???0.001/jmx1_events.fits'
    res1=findfile(name,count=n)
    res=[res,res1]
  endif 
  if (jemx eq 2) then begin 
    name='/home/sb/jemx/Crab/40_r3/00400002???0.001/jmx2_events.fits'
    res=findfile(name,count=n)
    name='/home/sb/jemx/Crab/40_r3/00400003???0.001/jmx2_events.fits'
    res1=findfile(name,count=n)
    res=[res,res1]
    name='/home/sb/jemx/Crab/40_r3/00400004???0.001/jmx2_events.fits'
    res1=findfile(name,count=n)
    res=[res,res1]
  endif

  
end

pro crab239files,res
  name='/home/sb/jemx/Crab/239/023900?????0.001/jmx?_events.fits'
  res=findfile(name,count=n)
  
end

pro crab422files,res
  name='/home/sb/jemx/Crab/422/042200?????0.001/jmx2_events.fits'
  res=findfile(name,count=n)
end

; get the time data
; INPUT   name   string containing the file name
;         jemx   JEM-X unit
; OUTPUT  time   array of time stamps in IJD of each detected photon 
pro getdata,name,jemx,time
  if (jemx eq 2) then fxbopen,unit,name,'JMX2-FULL-ALL'
  if (jemx eq 1) then fxbopen,unit,name,'JMX1-FULL-ALL'
  fxbread,unit,time,'TIME'
  fxbclose,unit
end

; get the time,x,y,pha data
pro getevdata,name,jemx,time,x,y,e
  if (jemx eq 2) then fxbopen,unit,name,'JMX2-FULL-ALL'
  if (jemx eq 1) then fxbopen,unit,name,'JMX1-FULL-ALL'
  fxbread,unit,time,'TIME'
  fxbread,unit,x,'RAWX'
  fxbread,unit,y,'RAWY'
  fxbread,unit,e,'PHA'
  fxbclose,unit
end

; read orbit information
; INPUT   o   orbit number
; OUTPUT  xyz   (3,N) array containing sat position in equatorial
; coordinates (km)
;        xyz    (3,N) array with sat velocity (km/s)
;        r     N array with sat distance from Earth center
;        v     N array with sat velocity (km/s)
pro readorbit,o,epoch,xyz,xyzvel,r,v,rdist,smaxis
  
  dname=string(o,format='(%"/r2/jemx/orbit/%4.4d/")')
  name=dname+'orbit_historic.fits'
  r= FILE_INFO(name)
  print,dname
  if (r.exists eq 0) then begin
     print,'Error find file: ',name
    return
  endif
  FXBOPEN, UNIT, name, 1
  FXBREAD, UNIT,epoch,'EPOCH'
  FXBREAD, UNIT,xyz,'XYZPOS'
  FXBREAD, UNIT,xyzvel,'XYZVEL'
  FXBREAD, UNIT,rdist,'RDIST'
  FXBREAD, UNIT,smaxis,'SMAXIS'
  FXBCLOSE, UNIT
  r=reform(sqrt(xyz(0,*)^2+xyz(1,*)^2+xyz(2,*)^2))
  v=reform(sqrt(xyzvel(0,*)^2+xyzvel(1,*)^2+xyzvel(2,*)^2))

end

; calculate sat projected velocity in particular direction
; INPUT   ra   R.A. of object in decimal degrees
;         dec  DEC of object in decimal degrees
;         xyzvel  (3,N) array of sat velocity in km/s
; OUTPUT  pv   N array of projected velocity in km/s
pro satv,ra,dec,xyzvel,pv 
; correct time to arrival at center of Earth
  ;ra=83.6332/!RADEG
  ;dec=22.0145/!RADEG
  rra=ra/!RADEG
  rdec=dec/!RADEG
  pv= xyzvel(0,*)*cos(rdec)*cos(rra) + $   ;Project velocity toward star
               xyzvel(1,*)*cos(rdec)*sin(rra) + xyzvel(2,*)*sin(rdec)
  
end

pro sat2e,tsat,xyz,s2e,ts,dt,res2  
; correct time to arrival at center of Earth
  ra=83.6332/!RADEG
  dec=22.0145/!RADEG
  s2e = xyz(0,*)*cos(dec)*cos(ra) + $   ;Project velocity toward star
               xyz(1,*)*cos(dec)*sin(ra) + xyz(2,*)*sin(dec)
  s2e=s2e/300000.0
  dt=INTERPOL(s2e,tsat,ts,/spline)
  ra=83.6332
  dec=22.0145
  res2= barycen(tsat+51544.0, ra, dec,orbit=xyz)
end

; conversion of sat time to barycenter arrival
; INPUT   xyz   (3,N) array giving sat position relative to Earth
; center
;         epoch  N array with the corresponding Earth time of xyz, IJD
;         time   array of the times stamps to be corrected, in IJD
; OUTPUT  epochc N array with the barycenter times of epoch
;         timebc array of barycenter corrected time stamps in IJD 
pro sat2barycenter,xyz,epoch,epochbc,time,timebc
  ; Crab ra and dec
  ra=83.6332D
  dec=22.0145D
  print,'test11'
  epochbc=barycen(epoch+51544.0,ra,dec,orbit=xyz)-51544.0
  print,'test2'
  timebc=INTERPOL(epochbc,epoch,time,/spline)
end

; convert IJD to a date (in TT reference)
; INPUT   ijd   value of IJD to be converted
; OUTPUT  yr,mn,day   year month and day of input
;         hr        hour of input as decimal number
pro ijd2date,ijd,yr,mn,day,hr
  offset=2451544.5D
  daycnv,ijd+offset,yr,mn,day,hr
  print,'IJD ',ijd,' is equal to: '
  print,yr,mn,day,hr,format='(i4.4,"  ",i2.2,"-",i2.2,"  " ,f9.6)'
  print,'yyyy  mo day  hour (decimal)'
end

; convert a date to IJD
pro date2ijd,yr,mn,day,hr,ijd
  offset=2451544.5D
  jdcnv,yr,mn,day,hr,jd
  ijd=jd-offset
end

; power density spectrum by FFT of event time stamps

pro powerd,t,nsample,f,pd
; INPUT
  ; t: time in seconds of each X-ray event
  ; nsample: number of samples per second of the light curve
; OUTPUT
;   f   frequency in Hertz
;   pd  power density of the FFT
;
  ; max power of 2
  p=long(2.0^findgen(25))
; generate a ligth curve fro the time stamps  
  h=histogram((t-t(0))*nsample)
  w=where(p lt n_elements(h))
  m=max(w)
  print,n_elements(h),m,p(m)
  h=h(0:p(m)-1)
  pd=abs(fft(h))
  pd=pd(0:p(m)/2-1)
  pd(0)=0
  ; normalize
  ;pd=pd*nsample*2.0/mean(h)
  print,'m=',mean(h)
  f=findgen(p(m)/2)/p(m)*nsample

end

; do period search on time stamps
; INPUT  
;  p0   center search period, in sec
;  dp   delta period for each step in search, s
;  n    number of steps in search around p0, to 2n+1
;  z    input time stamps of X-ray events in sec
; OUTPUT
;  p    2n+1 array of periods searched
;  f    merit function to measure deviation from flat signal
;  h    histogram of pulse from folding with the best period
;  pmax the best period derived from Gauss fit of f
;  dpmax the estimated sigma of the fitted pmax
;  a    the Gauss fit parameter array
;  sigma the sigma of a
;  r    2n+1 array of the gauss fit of f

pro psearch,p0,dp,n,z,p,f,h,pmax,dpmax,a,sigma,r
  f=dblarr(2*n+1)
  p=f
  nbins=20
  for i=-n,n do begin
    m=moment(histogram((z mod (p0+dp*i))/(p0+dp*i)*nbins))
    p(i+n)=p0+dp*i
    ;f(i+n)=sqrt(m(1))/m(0)*sqrt(n_elements(z)/20)
    ; variance divide by mean value over a trial priod
    ; if no periodicity at the trial period this ratio is 1
    f(i+n)=m(1)/m(0)
  endfor
  print,'max   period (s)'
  print,max(f,k),p(k)
  pmax=p(k)
  h=histogram((z mod (p0+dp*(k-n)))/(p0+dp*(k-n))*nbins)
  r=gaussfit(p,f,a,nterms=4,sigma=sigma)
  pmax=a(1)
  dpmax=sigma(1)
  
  print,'best period from gauss fit: ', pmax,' +-', a(2)

  plot,p,f,psym=10
  oplot,p,r,co=255
  h=histogram((z mod pmax)/(pmax)*100)
  
end

; print max of a function y of x in a range
; INPUT
; x   array of of X
; y    corresponding Y array
; x1  start of x range
; x2  end of x range
; OUTPUT
; val the x value where y has maximum in the x1,x2 range
pro printmax,x,y,x1,x2,val
  w=where(x gt x1 and x lt x2,n)
  if (n lt 1) then print,"illegal range"
  if (n gt 1) then print,max(double(y(w)),n),double(x(w(n))),format='("max(y) = ",e15.8,"  for x = ",e15.8)'  
  val=double(x(w(n)))
end

; plot a pulse pulse shape repeated over two periods
; INPUT  h  array with pulse histogram over one period
; OUTPUT phase  phase of output pulse, range 0,2
; hh   the pulse repeated
;
pro plotpulse,h,phase,hh
  n=n_elements(h)
  x=findgen(2*n)/n
  print,max(h,m),m
  hh=shift([h,h],-m+n/2)
  w=where(x gt 0.1 and x lt 0.35)
  
  y0=mean(hh(w))
  plot,x,hh/y0,psym=10,ysty=1,tit='pulse profile (plotted twice) normalized to 1 at minimum flux',xtit='phase'  
  
end


; save the IDL plot window as a .png graphics file
; INPUT  name  string containing the output graphics file name
; 
; (note than the extension .png should be given in the string)
; the graphics color table is inverted to produce black on white plots
; from the scree plot of white on black
;
pro savepng,name
  im=tvrd()
  write_png,name,255B-im
end 

; do up to all included in step 5
; return the time in t and the barycentered time in tbc
pro step5,t,tbc
  step2,name,time,t,x,y,e,im
  readorbit,40,epoch,xyz,xyzvel,r,v
  sat2barycenter,xyz,epoch,epochbc,time,timebc
  tbc=(timebc-timebc(0))*86400.0
end