! @(#)airmass.prg	17.1.1.1 (ESO-DMD) 01/25/02 17:12:18
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.NAME  airmass.prg
!
!.AUTHOR  M. Rosa (ST_ECF)	VERSION 2.0	29-MAR-1986
!         D. Baade (ST-ECF)  870207
!	  KB  901108, 920226, 920319, 941128
!         PJG 920904
!
!.KEYWORDS
! airmass calculation, Julian date, conversion between UT and ST
!							
!.PURPOSE  Calculate airmass ( = f(sec z) ) and JD and sidereal time
!          Use and update descriptors in frame option.
!          Order of positional parameters longitude and latitude exchanged
!          convert to eastern longitude!  PJG
!          Input has been changed so that eastern longitude is required
!          following IAU standards.     PJG 
!
!.CALL  a) AIRMASS image [long] [lat]
!       b) AIRMASS ra dec [ST] [dur] [long] [lat] date UT
!
! a) input from and to image descriptors
!
!           The descriptor O_TIME is organized as follows:
!
!           1) YEAR  2) MONTH  3) DAY  4) JD  5) UT (hrs...)
!           6) ST (hrs...)  7) DUR (seconds) 8) Offset of UT clock (hrs...)
!
!           ST and JD will always be (re-)calculated by the program.
!           (Use command WRITE/DESC if UT, date, and/or duration are not set)
!
!           O_POS(1): right ascension in real degrees (!!, not: hours)
!           O_POS(2): declination in real degrees
!
! b) input only from keywords
!
!	    ra     right ascension in HH,MM,SS
!	    dec    declination in DD,MM,SS
!	    ST     sidereal time in HH,MM,SS (DEFAULT: = RA)
!                  at half exposure
!                  if ST > 25 ---> ST = ra
!                  if ST < 0, ST will be calculated from UT and date
!           dur    exposure time (default: 0.0)
!           long   east longitude in DD,MM,SS.ss (DEFAULT: ESO -70,43,55.35)
!	    lat    latitude in DD,MM,SS.ss  (DEFAULT: ESO -29,15,25.8)
!           date   year,month day
!           UT     UT in HH,MM,SS.ss..
!
!           If DUR > 0.0, some mean airmass for exposure time is calculated
!             (arithmetic mean of AIRM(ST), AIRM(ST+0.5DUR), AIRM(ST+DUR))
!
!------------------------------------------------------------------------------
!
!.INPUT (interactive mode)
!      Right ascension (HH,MM,SS.ss...) INPUTR(1:3)
!      Declination (DD,MM,SS.ss...)     INPUTR(4:6)
!      Sidereal time (HH,MM,SS.ss...)   INPUTR(7:9)
!      exposure time (sec)              INPUTR(10:10)
!      East longitude (DD,MM,SS.ss...)  INPUTR(11:13)
!      Latitude (DD,MM,SS.ss...)        INPUTR(14:16)
!      UT (HH,MM,SS.ss...)              INPUTR(17:19)
!      year,month,day                   INPUTI(1:3)
!
!      If exposure time > 0.0 sec the airmass calculated is the
!      average of airmasses at t0, 0.5(t1+t0) and t1
!
!      UT+longitude+year,month,day and ST are equivalent.
!      Therefore: non-positive SID resp. non-positive UT request the
!                 UT resp. SID to be used in calc. of airmass
!
!.OUTPUT  Airmass
!
!---------------------------------------------------------------
!
define/local signs/c/1/2 ++	           ! signs of latitude and declination
define/local rneg/r/1/1 1.0
!
set/format g10.5,g24.12
!
define/param p1 ? ? "Enter image or right ascension (hour,min,sec): "
!
if m$index(p1,",") .ne. 0 goto keys	!check, if frame or number(s)
if m$tstno(p1) .ne. 0 goto keys
!
! ... a) input and output from and to descriptors of image
!
define/param p1 + ima 
!
define/local image/c/1/100 {p1}		!name of image
define/local long/c/1/40 0		!eastern longitude
define/local lat/c/1/40  0		!latitude
!  
define/param p2 0 N  ! eastern longitude of observatory (no input=0)
define/param p3 0 N  ! latitude of observatory          (no input=0)
!
! check first if input for longitude/latitude on line:
!
if  p2 .eq. "0" .and. p3 .eq. "0" then
  ! no input on command line found.
  ! check if descriptors there
  if m$existd(p1,"ESO.TEL.GEOLON") .eq. 1 .AND. m$existd(p1,"ESO.TEL.GEOLAT") .eq. 1 then
     write/out ** Descriptors for longitude and latitude of telescope found **
     write/out ** ESO.TEL.GEOLON and ESO.TEL.GEOLAT will be used **
     ! 
     write/key long/c/1/40 {{p1},ESO.TEL.GEOLON},0,0
     write/key lat/c/1/40 {{p1},ESO.TEL.GEOLAT},0,0
  else
     write/out ** No descriptors ESO.TEL.GEOLON, ESO.TEL.GEOLAT found **
     ! if no long and lat AND no descriptors there, use default: La Silla
     p2 = "-70.,43.,55.35"
     p3 = "-29.,15.,25.8"
     write/out ** Default (=longitude and latitude of La Silla) will be used **
     write/key long/c/1/40 {p2},0,0	!eastern longitude
     write/key lat/c/1/40  {p3},0,0	!latitude
  endif
else
  ! yes, input on command line:
  write/out ** Longitude: {p2} and latitude: {p3} used **  
  write/key long/c/1/40 {p2},0,0	!eastern longitude
  write/key lat/c/1/40  {p3},0,0	!latitude
endif
!
! ... first compute ST and JD
! ... (results in keywords OUTPUTD(19) and OUTPUTD(20) as well as
! ... image descriptors O_TIME(6) and O_TIME(4) respectively
!
compute/st {image} {long}
!
inputr(1) = outputd(19)*15.-m$value({image},o_pos(1)) !H.A. = ST-R.A. (hour_ang)
write/keyw inputr/r/2/3 {{image},o_pos(2)},0,0	 !declination
if inputr(2) .lt. 0.0 signs(2:2) = "-"	         !negative declination	
inputr(5) = m$value({image},o_time(7))           !length of exposure
write/keyw inputr/r/6/3 {lat}		         !latitude of observatory
signs(1:1) = lat(1:1)
rneg = (inputr(6)*3600.)+(inputr(7)*60.)+inputr(8)
if rneg .lt. 0.0 signs(1:1) = "-"
!
run MID_EXE:AIRMASS			       ! compute mean airmass
!
copy/kd outputr/r/1/1 {p1} o_airm/r/1/1        ! write result to descriptor
!
write/out "Image:    {IMAGE}          Mean airmass:      {OUTPUTR(1)}"
!
return					       ! job done
!
KEYS:
!
! Interactive option:
!
define/param p2 ? ? "Enter declination (degree,min,sec):         "
define/param p3 -1.0. N			       ! ST
define/param p4 0.0 N                          ! exposure time
define/param p5 -70.,43.,55.35 N               ! east longitude (def. La Silla)
define/param p6 -29.,15.,25.8 N                ! latitude (def.: La Silla)
!
define/local alpha0/r/1/3 {p1},0,0             ! right ascension as hour,min,sec
define/local alpha/r/1/1 0                     ! right ascension in hours
define/local delta/c/1/40 "{p2},0,0"           ! declination
define/local st/r/1/1 {p3}		       ! ST
define/local exptim/r/1/1 {p4}                 ! exposure time (in seconds)
define/local long/c/1/40 "{p5},0,0"	       ! longitude
define/local lat/c/1/40 "{p6},0,0"             ! latitude
!
! ... if ST not known (i.e. entered as negative) but wanted calculate it first
!
if st .lt. 0 then
   define/param p7 ? N "Enter date of observation (year,month,day): "
   define/param p8 ? N "Enter UT of observation (hour,min,sec):     "
!
   define/local date/c/1/40 "{p7},0,0"	     ! date
   define/local ut/c/1/40 "{p8},0,0"	     ! UT
!
   compute/st {date} {ut} {long}
   st = outputd(19)
endif
!
alpha = alpha0(1)+(alpha0(2)+alpha0(3)/60.)/60.    ! R.A. in real hours
!
! ... other special option: only calculate minimal airmass,
! ... i.e. for object on meridian
!
if st .gt. 25. st = alpha
!
inputr(1) = (st-alpha)*15.                       ! H.A. = ST - R.A. (hour angle
write/keyw inputr/r/2/3 {delta},0,0	         ! declination
rneg = (inputr(2)*3600.)+(inputr(3)*60.)+inputr(4)
signs(2:2) = delta(1:1)
if rneg .lt. 0.0 signs(2:2) = "-"		 ! negative declination
inputr(5) = exptim		   	         ! length of exposure
write/keyw inputr/r/6/3 {lat},0,0                 ! latitude of observatory
signs(1:1) = lat(1:1)
rneg = (inputr(6)*3600.)+(inputr(7)*60.)+inputr(8)
if rneg .lt. 0.0 signs(1:1) = "-" 		 ! negative latitude
!
run MID_EXE:AIRMASS				         ! compute airmass
!
write/out  "Airmass:                     {OUTPUTR(1)}"	     ! display result