; docformat = 'rst'
;+
; :Description:
; Calculates rise, transit and set times, plus the time above the horizon (sky time) for the given
; equatorial coordinates, location and date.
;
; GDL-compatible. All dependencies are routines from idlastro (http://idlastro.gsfc.nasa.gov/).
;
; :Returns:
; A structure with the fields transit, sky, rise and set, each with the same dimensions as ra,dec, with
; the times for transit, above the horizon, rise and set, respectively, in civil hours. Those points that are never
; above the horizon have sky=0 and NaN for rise and set. Those points that are always above the horizon have sky=24
; and NaN for rise and set.
;
; :Params:
; ra : in, required
; Right ascension(s) of the point(s) to calculate the visibility for. Units: hours. May be a scalar or arrar,
; and must have the same number of elements as dec.
; dec : in, required
; Declination(s) of the point(s) to calculate the visibility for. Units: degrees. May be a scalar or arrar,
; and must have the same number of elements as ra.
;
; :Keywords:
; obsname : in, optional
; String with the observatory name, as understood by idlastro's observatory procedure. If provided, supersedes lat, lon and tz.
; lat : in, optional
; Latitude of the observing location, in degrees (positive is north).
; lon : in, optional
; Longitude of the observing location, in degrees (positive is east, contrary to that returned by observatory).
; tz : in, optional
; Time zone of the observing location, in hours (positive is east, contrary to that returned by observatory).
; day : in, required
; Day of the month for the observation.
; month : in, required
; Month for the observation.
; year : in, required
; Year for the observation.
; prev_day : in, optional, default=0
; By default, returned times of day (for transit, rise and set) are relative to local midnight, so times before
; midnight are negative. If this keyword is set, times before midnight are returned as the time for the previous day,
; instead.
; extra : in, optional, default=0
; If set, extra data is output (sunrise/sunset times, target altitude at sunrise/sunet, etc).
;
; :Examples:
;
; Make up some coordinates and calculate the times::
;
; ra=[22d0+3d0/60d0+18d0/3600d0,0d0,19d0]
; dec=[-9d0-40d0/60d0-34d0/3600d0,-87d0,80d0]
; r=pp_rise_set(ra,dec,obsname='eso',day=0,month=06,year=2010,/prev)
; print,r.transit
; 6.2173759 8.1019913 3.1540115
; print,r.sky
; 12.765635 24.000000 0.0000000
; print,r.rise
; 23.834559 NaN NaN
; print,r.set
; 12.600193 NaN NaN
;
; :Uses: observatory, ten, ct2lst, jdcnv
;
; :Author: Paulo Penteado (pp.penteado@gmail.com), Jun/2010
;-
function pp_rise_set,ra,dec,obsname=obsname,lat=lat,lon=lon,tz=tz,day=day,$
month=month,year=year,prev_day=prev,extra=extra
compile_opt idl2,logical_predicate
;Defaults
extra=keyword_set(extra)
prev=n_elements(prev) eq 1 ? prev : 0
if (n_elements(day) ne 1) || (n_elements(month) ne 1) || (n_elements(year) ne 1) then $
message,'Time must be provided by day,month,year'
if (n_elements(obsname) eq 1) then begin
observatory,obsname,obs
lat=obs.latitude
lon=-obs.longitude
tz=-obs.tz
endif else if (n_elements(lat) ne 1) || (n_elements(lon) ne 1) || (n_elements(tz) ne 1) then $
message,'Location must be provided by either obsname, or lat,lon,tz'
;Find the transit time
;Find the LST at 0h for the day
ct2lst,lst,lon,tz,0d0,day,month,year
;Find the LST at 0h for the previous day
ct2lst,lst1,lon,tz,0d0,day-1,month,year
;Sidereal to Civil time factor
scfac=1d0+(lst-lst1)/24d0
;Transit time in sidereal hours after midnight
transit_sid=ra-lst
;Transit time in civil hours after midnight
transit_civ=transit_sid*scfac
;Rise/set/sky time
rise=double(ra) & rise[*]=0d0 & set=rise & sky=rise
;Those that are never above the horizon
w=lat gt 0d0 ? where(dec lt (-90d0+lat),nw) : where(dec gt 90d0+lat,nw)
if (nw gt 0) then begin
rise[w]=!values.d_nan & set[w]=!values.d_nan & sky[w]=0d0
endif
;Those that are always above the horizon
w=lat gt 0d0 ? where(dec gt (90d0-lat),nw) : where(dec lt -90d0-lat,nw)
if (nw gt 0) then begin
rise[w]=!values.d_nan & set[w]=!values.d_nan & sky[w]=24d0
endif
;Those that rise and set
w=where(finite(rise),nw)
if (nw gt 0) then begin
;Find the hour angle of set time
rdec=dec[w]*!dpi/180d0 & rlat=lat[w]*!dpi/180d0
tmp=dcomplex(-sin(rlat)*sin(rdec)/(cos(rlat)*cos(rdec)))
ha=acos(real_part(tmp))*12d0/!dpi
;Time in the sky (civil hours)
sky[w]=2d0*ha*scfac
;Rise time (civil hours after midnight)
rise[w]=transit_civ-ha*scfac
;Set time (civil hours after midnight)
set[w]=transit_civ+ha*scfac
endif
;Convert negative times to previous day
if (prev) then begin
rise=(rise+24d0) mod 24
set=(set+24d0) mod 24
transit_civ=(transit_civ+24d0) mod 24
endif
if extra then begin
jd=julday(month,day,year)
sunpos,jd,sunra,sundec
sunra/=15d0
sun=pp_rise_set(sunra,sundec,obsname=obsname,lat=lat,lon=lon,tz=tz,day=day,$
month=month,year=year,prev_day=prev)
sun=create_struct(sun,'ra',sunra,'dec',sundec)
eq2hor,ra*15d0,dec,jd+(sun.rise+tz+12d0)/24d0,altr,azr,lat=lat,lon=lon
eq2hor,ra*15d0,dec,jd+(sun.set+tz+12d0)/24d0,alts,azs,lat=lat,lon=lon
target_nosun={sunset_alt:alts,sunrise_alt:altr};,no_sun_hours:hours}
endif
;Pack the results
ret={transit:transit_civ,sky:sky,rise:rise,set:set}
if extra then begin
if n_elements(obsname) eq 0 then obsname=''
ret=create_struct(ret,'jd',jd,'sun',sun,'target',target_nosun,'month',$
month,'day',day,'year',year,'ra',ra,'dec',dec,'obsname',obsname,'lat',lat,'lon',lon,'tz',tz)
endif
return,ret
end