; -------------------------------------------------------------------- FUNCTION nlines,file ; This function opens the specified file, determines its length in ; bytes, and reads it in as a BYTE-valued array. It then returns ; the number of occurrences of 10B in the file. 10B is the ASCII ; code for a line feed, which marks the end of a line of text (in ; UNIX, at least). ; Modified 10 January 1995 by D P Steele: ; If called under VMS, this function dodges the whole issue of ; counting lines and SPAWNs a call to WC.COM, which uses EVE (the ; Extensible Vax Editor) to count the number of lines and return ; the result. Obviously, WC.COM must be accessible to the function. IF !VERSION.OS EQ 'vms' THEN BEGIN SPAWN,'@$DISC1:[STEELE.IDL]WC '+file,n n=LONG(n(0)) ENDIF ELSE BEGIN OPENR,unit,file,/GET_LUN f=FSTAT(unit) IF f.SIZE EQ 0 THEN n=0L ELSE BEGIN b=BYTARR(f.size) READU,unit,b n=LONG(TOTAL(b EQ 10B)) ENDELSE CLOSE,unit FREE_LUN,unit ENDELSE RETURN,n END ;------------------------------------------------------------- ;+ ; NAME: ; QGETFILE ; PURPOSE: ; Read a text file into a string array. ; CATEGORY: ; CALLING SEQUENCE: ; s = qgetfile(f) ; INPUTS: ; f = text file name. in ; KEYWORD PARAMETERS: ; Keywords: ; ERROR=err error flag: 0=ok, 1=file not opened. ; /QUIET means give no error message. ; OUTPUTS: ; s = string array. out ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner, 20 Mar, 1990 ; D. Steele, 18 Mar, 1994 - added line count, changed WHILE loop ; to FOR loop, and changed name from ; GETFILE.PRO to QGETFILE.PRO. ; ; Copyright (C) 1990, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- function qgetfile, file, error=err, help=hlp, quiet=quiet if (n_params(0) lt 1) or keyword_set(hlp) then begin print,' Read a text file into a string array.' print,' s = getfile(f)' print,' f = text file name. in' print,' s = string array. out' print,' Keywords:' print,' ERROR=err error flag: 0=ok, 1=file not opened.' print,' /QUIET means give no error message.' return, -1 endif get_lun, lun on_ioerror, err ; The following lines were added by D P Steele, 94/03/18 nl=nlines(file) s=STRARR(nl) openr, lun, file ; s = [' '] t = '' ; while not eof(lun) do begin FOR i=0L,nl-1 DO BEGIN ; replace EOF search with FOR loop readf, lun, t ; s = [s,t] s(i)=t ; avoid slow string concatenation ; endwhile ENDFOR close, lun free_lun, lun err = 0 ; return, s(1:*) RETURN,s err: if !err eq -168 then begin if not keyword_set(quiet) then print,' Non-standard text file format.' free_lun, lun ; return, s(1:*) RETURN,s endif if not keyword_set(quiet) then print,$ ' Error in getfile: File '+file+' not opened.' free_lun, lun err = 1 return, -1 end ; -------------------------------------------------------------------- ;+ ; NAME: ; MPUTFILE ; PURPOSE: ; Write a text file from a string array. ; CATEGORY: ; CALLING SEQUENCE: ; mputfile, f, s ; INPUTS: ; f = text file name. in ; s = string array. in ; KEYWORD PARAMETERS: ; Keywords: ; ERROR=err error flag: 0=ok, 1=invalid string array. ; OUTPUTS: ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner, 20 Mar, 1990 ; R. Sterner, 4 Nov, 1992 --- allowed scalar strings. ; ; Copyright (C) 1990, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- pro mputfile, file, s, error=err, help=hlp if (n_params(0) lt 1) or keyword_set(hlp) then begin print,' Write a text file from a string array.' print,' mputfile, f, s' print,' f = text file name. in' print,' s = string array. in' print,' Keywords:' print,' ERROR=err error flag: 0=ok, 1=invalid string array.' return endif if datatype(s) ne 'STR' then begin print,' Error in mputfile: argument must be a string array.' err = 1 return endif get_lun, lun openw, lun, file for i = 0L, n_elements(s)-1 do begin t = s(i) if t eq '' then t = ' ' printf, lun, t endfor close, lun free_lun, lun err = 0 return end ; -------------------------------------------------------------------- ;+ ; NAME: ; MONTHNAMES ; PURPOSE: ; Returns a string array of month names. ; CATEGORY: ; CALLING SEQUENCE: ; mnam = monthnames() ; INPUTS: ; KEYWORD PARAMETERS: ; OUTPUTS: ; mnam = string array of 13 items: out ; ['Error','January',...'December'] ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner, 18 Sep, 1989 ; ; Copyright (C) 1989, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- function monthnames, help=hlp if keyword_set(hlp) then begin print,' Returns a string array of month names.' print,' mnam = monthnames()' print,' mnam = string array of 13 items: out' print," ['Error','January',...'December']" return, -1 endif mn = ['Error','January','February','March','April','May',$ 'June','July','August','September','October',$ 'November','December'] return, mn end ; -------------------------------------------------------------------- ;+ ; NAME: ; YMD2JD ; PURPOSE: ; From Year, Month, and Day compute Julian Day number. ; CATEGORY: ; CALLING SEQUENCE: ; jd = ymd2jd(y,m,d) ; INPUTS: ; y = Year (like 1987). in ; m = month (like 7 for July). in ; d = month day (like 23). in ; KEYWORD PARAMETERS: ; OUTPUTS: ; jd = Julian Day number (like 2447000). out ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner, 23 June, 1985 --- converted from FORTRAN. ; Johns Hopkins University Applied Physics Laboratory. ; RES 18 Sep, 1989 --- converted to SUN ; ; Copyright (C) 1985, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- function ymd2jd, iy, im, id, help=hlp if (n_params(0) LT 3) or keyword_set(hlp) then begin print,' From Year, Month, and Day compute Julian Day number.' print,' jd = ymd2jd(y,m,d)' print,' y = Year (like 1987). in' print,' m = month (like 7 for July). in' print,' d = month day (like 23). in' print,' jd = Julian Day number (like 2447000). out' return, -1 endif y = long(iy) m = long(im) d = long(id) jd = 367*y-7*(y+(m+9)/12)/4-3*((y+(m-9)/7)/100+1)/4 $ +275*m/9+d+1721029 return, jd end ; -------------------------------------------------------------------- ;+ ; NAME: ; WEEKDAY ; PURPOSE: ; Compute weekday given year, month, day. ; CATEGORY: ; CALLING SEQUENCE: ; wd = weekday(y,m,d,[nwd]) ; INPUTS: ; y, m, d = Year, month, day (Like 1988, 10, 31). in ; KEYWORD PARAMETERS: ; OUTPUTS: ; wd = Returned name of weekday. out ; nwd = optional Weekday number. out ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner. 31 Oct, 1988. ; Johns Hopkins University Applied Physics Laboratory. ; RES 18 Sep, 1989 --- converted to SUN ; ; Copyright (C) 1988, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- FUNCTION WEEKDAY, Y, M, D, NWD, help=hlp IF (N_PARAMS(0) LT 3) or keyword_set(hlp) THEN BEGIN PRINT,' Compute weekday given year, month, day.' PRINT,' wd = weekday(y,m,d,[nwd])' PRINT,' y, m, d = Year, month, day (Like 1988, 10, 31). in' PRINT,' wd = Returned name of weekday. out' PRINT,' nwd = optional Weekday number. out' RETURN, -1 ENDIF NAMES = ['','Sunday','Monday','Tuesday',$ 'Wednesday','Thursday','Friday','Saturday'] NWD = ((YMD2JD(Y,M,D) + 1) MOD 7) + 1 RETURN, NAMES(NWD) END ; -------------------------------------------------------------------- ;+ ; NAME: ; STREP ; PURPOSE: ; Edit a string by position. Precede, Follow, Replace, Delete. ; CATEGORY: ; CALLING SEQUENCE: ; newstring = strep(string,cmd,p,ss,[iflg]) ; INPUTS: ; string = string to edit. in ; cmd = edit command: in ; 'P' = precede position p with substring ss. ; 'F' = follow position p with substring ss. ; 'R' = replace text starting at position p ; with text from substring ss. ; 'D' = delete N characters starting at ; position p. The calling sequence for ; this command is slightly different: ; IFLG = STREP(string,'D',p,n,[iflg]) ; Where n = number of characters to delete. ; p = character position to use. in ; 0 = first char. Any number larger ; than the string length = last char. ; ss = substring to use. For 'D' command in ; n is used instead of ss. ; KEYWORD PARAMETERS: ; OUTPUTS: ; iflg = 0 for a successful edit, out ; iflg = -1 for an error and no change to string. ; newstring = edited string. out ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; Written by R. Sterner, 27 Dec, 1984. ; Converted to SUN 13 Aug, 1989 --- R. Sterner. ; Johns Hopkins University Applied Physics Laboratory. ; ; Copyright (C) 1984, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- FUNCTION strep,S,CMD,IP,SS,IFLG, help = h IF (N_PARAMS(0) LT 4) or keyword_set(h) THEN BEGIN PRINT,' Edit a string by position. Precede, Follow, Replace, Delete.' print,' newstring = strep(string,cmd,p,ss,[iflg]) print,' string = string to edit. in' print,' cmd = edit command: in' print," 'P' = precede position p with substring ss. print," 'F' = follow position p with substring ss. print," 'R' = replace text starting at position p print,' with text from substring ss. print," 'D' = delete N characters starting at print,' position p. The calling sequence for print,' this command is slightly different: print," IFLG = STREP(string,'D',p,n,[iflg]) print,' Where n = number of characters to delete. print,' p = character position to use. in' print,' 0 = first char. Any number larger print,' than the string length = last char. print," ss = substring to use. For 'D' command in" print,' n is used instead of ss. print,' iflg = 0 for a successful edit, out' print,' iflg = -1 for an error and no change to string. print,' newstring = edited string. out' RETURN, -1 ENDIF N = STRLEN(S) - 1 IF N LT 0 THEN $ IF CMD EQ 'D' THEN RETURN, '' ELSE RETURN, SS P = IP>00) EQ SSL-3) THEN SSTYP = 2 IF (PDOT GT 0) AND (PDOT EQ SSL-3) THEN SSTYP = 2 IF (PDOT EQ 0) AND (SSL GT 3) THEN SSTYP = 3 IF (PDOT GT 0) AND (PDOT LT SSL-3) THEN SSTYP = 4 IF (PDOT EQ 0) AND (SSL EQ 3) THEN SSTYP = 5 NED = 0 ; number of occurrences actually changed. CASE SSTYP OF 1: BEGIN S = OLD E = '' END 2: BEGIN S = STRSUB(OLD,0,SSL-4) E = '' END 3: BEGIN S = '' E = STRSUB(OLD,3,SSL-1) END 4: BEGIN S = STRSUB(OLD,0,PDOT-1) E = STRSUB(OLD,PDOT+3,SSL-1) END 5: BEGIN S = '' E = '' END ELSE: PRINT, 'ERROR IN SSTYP' ENDCASE ;--------------- Find substring # N start --------------- POS = -1 nfor = n>1 LOOP: FOR I = 1, nfor DO BEGIN POS = POS + 1 CASE SSTYP OF 1: POS = STRPOS(RSTR,S,POS) 2: POS = STRPOS(RSTR,S,POS) 4: POS = STRPOS(RSTR,S,POS) 3: POS = STRPOS(RSTR,E,POS) 5: POS = 0 ENDCASE IF POS LT 0 THEN GOTO, DONE ENDFOR ;---------- Find substring # N END ---------------- CASE SSTYP OF 1: BEGIN POS1 = POS POS2 = POS + STRLEN(S) - 1 END 2: BEGIN POS1 = POS POS2 = STRLEN(RSTR) - 1 END 3: BEGIN POS1 = 0 POS2 = POS + STRLEN(E) - 1 END 4: BEGIN POS1 = POS POS2 = STRPOS(RSTR,E,POS+1) IF (POS2 LT 0) THEN GOTO, DONE POS2 = POS2 + STRLEN(E) - 1 END 5: BEGIN POS1 = 0 POS2 = STRLEN(RSTR) - 1 END ENDCASE ;------------ edit string -------------- CASE CMD OF 'P': BEGIN RSTR = STREP(RSTR,CMD,POS1,NEW) POS = POS + STRLEN(NEW) END 'F': BEGIN RSTR = STREP(RSTR,CMD,POS2,NEW) POS = POS + STRLEN(NEW) END 'R': BEGIN RSTR = STREP(RSTR,'D',POS1,POS2-POS1+1) IF (POS1 GT 0) THEN $ RSTR = STREP(RSTR,'F',POS1-1,NEW) IF (POS1 EQ 0) THEN $ RSTR = STREP(RSTR,'P',0,NEW) POS = POS + STRLEN(NEW) - 1 END 'D': BEGIN RSTR = STREP(RSTR,CMD,POS1,POS2-POS1+1) POS = POS - 1 END ELSE: BEGIN PRINT, 'Error in cmd' RETURN,RSTR END ENDCASE NED = NED + 1 IF SSTYP EQ 5 THEN RETURN,RSTR IF N EQ 0 THEN GOTO, LOOP DONE: RETURN, RSTR END ; -------------------------------------------------------------------- ;+ ; NAME: ; YMD2DATE ; PURPOSE: ; Convert from year, month, day numbers to date string. ; CATEGORY: ; CALLING SEQUENCE: ; date = ymd2date(Y,M,D) ; INPUTS: ; y = year number (like 1986). in ; m = month number (1 - 12). in ; d = day of month number (1 - 31). in ; KEYWORD PARAMETERS: ; Keywords: ; FORMAT = format string. Allows output date to be customized. ; The following substitutions take place in the format string: ; Y$ = 4 digit year. ; y$ = 2 digit year. ; N$ = full month name. ; n$ = 3 letter month name. ; d$ = day of month number. ; W$ = full weekday name. ; w$ = 3 letter week day name. ; OUTPUTS: ; date = returned date string (like 24-May-1986). out ; COMMON BLOCKS: ; NOTES: ; Notes: ; The default format string is 'd$-n$-Y$' giving 24-Sep-1989 ; Example: FORMAT='w$ N$ d$, Y$' would give 'Mon ; MODIFICATION HISTORY: ; R. Sterner. 16 Jul, 1986. ; RES 18 Sep, 1989 --- converted to SUN ; R. Sterner, 28 Feb, 1991 --- modified format. ; R. Sterner, 16 Dec, 1991 --- added space to 1 digit day. ; Johns Hopkins University Applied Physics Laboratory. ; ; Copyright (C) 1986, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- function ymd2date, y, m, d, help=hlp, format=frmt if (n_params(0) lt 3) or keyword_set(hlp) then begin print,' Convert from year, month, day numbers to date string.' print,' date = ymd2date(Y,M,D)' print,' y = year number (like 1986). in' print,' m = month number (1 - 12). in' print,' d = day of month number (1 - 31). in' print,' date = returned date string (like 24-May-1986). out' print,' Keywords:' print,' FORMAT = format string. Allows output date to be '+$ 'customized.' print,' The following substitutions take place in the '+$ 'format string:' print,' Y$ = 4 digit year.' print,' y$ = 2 digit year.' print,' N$ = full month name.' print,' n$ = 3 letter month name.' print,' d$ = day of month number.' print,' W$ = full weekday name.' print,' w$ = 3 letter week day name.' print,' Notes:' print," The default format string is 'd$-n$-Y$' giving 24-Sep-1989" print," Example: FORMAT='w$ N$ d$, Y$' would give 'Mon "+$ "September 18, 1989'" return, -1 endif ;---- error check ----- w = where(y lt 0, cnt) if cnt gt 0 then begin print,'Error in ymd2date: invalid year: ',y(w) return, -1 endif ;----- Handle 2 digit years ------ w = where(y lt 50, cnt) ; Y < 50 assumed 2xxx. if cnt gt 0 then y(w) = y(w) + 2000 w = where(y lt 100, cnt) ; 50 < Y < 100 assumed 19xx. if cnt gt 0 then y(w) = y(w) + 1900 w = where((m lt 1) or (m gt 12),cnt) if cnt gt 0 then begin print,'Error in ymd2date: invalid month: ',m(w) return, -1 endif w = where((d lt 1) or (d gt 31), cnt) if cnt gt 0 then begin print,'Error in ymd2date: invalid month day: ', d(w) return, -1 endif ;----- format string ------ fmt = 'd$-n$-Y$' if keyword_set(frmt) then fmt = frmt n_dates = n_elements(y) date_arr = strarr(n_dates) ;----- Loop through all dates --------- for i = 0, n_dates-1 do begin ;----- Get all the allowed parts ----- yu = strtrim(y(i),2) yl = strtrim(fix(y(i)-100*fix(y(i)/100)),2) mnames = monthnames() mu = mnames(m(i)) ml = strmid(mu,0,3) dl = strtrim(d(i),2) if strlen(dl) eq 1 then dl = ' '+dl ; Add space to 1 digit day. wu = weekday(y(i),m(i),d(i)) wl = strmid(wu,0,3) ;---- Do substitutions ------ date = fmt date = stress(date, 'R', 0, 'Y$', yu) ; 4 digit year. date = stress(date, 'R', 0, 'y$', yl) ; 2 digit year. date = stress(date, 'R', 0, 'd$', dl) ; day of month. date = stress(date, 'R', 0, 'N$', mu) ; Full month name. date = stress(date, 'R', 0, 'n$', ml) ; 3 letter month name. date = stress(date, 'R', 0, 'W$', wu) ; Full weekday name. date = stress(date, 'R', 0, 'w$', wl) ; 3 letter weekday name. date_arr(i) = date endfor if n_dates eq 1 then date_arr=date_arr(0) ; Return scalars as scalars. return, date_arr end ; -------------------------------------------------------------------- ;+ ; NAME: ; SETMINUS ; PURPOSE: ; Eliminate elements from a set that are also in another set. ; CATEGORY: ; CALLING SEQUENCE: ; c = setminus(a,b) ; INPUTS: ; a = set (array) to process. in ; b = set (array) of elements to eliminate. in ; KEYWORD PARAMETERS: ; OUTPUTS: ; c = resulting set (array). out ; COMMON BLOCKS: ; NOTES: ; Notes: works for any type arrays. ; Will not remove the element of a one element array. ; This is because IDL does not have empty arrays ; (corresponding to empty sets). ; MODIFICATION HISTORY: ; R. Sterner, 7 Aug, 1991 ; ; Copyright (C) 1991, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- function setminus, a, b, help=hlp if (n_params(0) lt 2) or keyword_set(hlp) then begin print,' Eliminate elements from a set that are also in another set.' print,' c = setminus(a,b)' print,' a = set (array) to process. in' print,' b = set (array) of elements to eliminate. in' print,' c = resulting set (array). out' print,' Notes: works for any type arrays.' print,' Will not remove the element of a one element array.' print,' This is because IDL does not have empty arrays' print,' (corresponding to empty sets).' return, -1 endif c = a for i = 0, n_elements(b)-1 do begin w = where(c ne b(i), count) if count gt 0 then c = c(w) endfor return, c end ; -------------------------------------------------------------------- FUNCTION PoCasite,year,month,day yr=year MOD 100 ns=N_ELEMENTS(yr) PCsite=REPLICATE(1,ns) ; default location is Calgary ; Decide where Polar Camera was located rabb=WHERE((yr EQ 93) AND ((month EQ 1) OR (month EQ 2)),nrabb) IF nrabb GT 0 THEN PCsite(rabb)=2 ; Rabbit Lake in Jan, Feb 1993 eur=WHERE(((yr EQ 93) AND (month GE 10)) OR (yr GE 94),neur) IF neur GT 0 THEN PCsite(eur)=3 ; Eureka since October 1993 IF ns EQ 1 THEN PCsite=PCsite(0) ; scalar inputs so scalar output RETURN,PCsite END ; -------------------------------------------------------------------- FUNCTION ymdh2jdf,year,month,day,hour ; This function is basically a wrapper for ymd2jd.pro from the ; JHU/APL IDL library. It includes the specification of a (possibly ; fractional) hour of the day, and returns a (possibly fractional) ; Julian date. ijd=ymd2jd(year,month,day) RETURN,DOUBLE(ijd)+(hour-12.)/24. END ; -------------------------------------------------------------------- PRO MOONPOS, dte, ra, dec, x, y, z ;+ ; NAME: ; MOONPOS ; PURPOSE: ; To compute the RA and Dec of the Moon at a given date. ; ; CALLING SEQUENCE: ; MOONPOS, Dte, Ra, Dec, [X, Y, Z] ; ; INPUTS: ; Dte - The Julian date of the day (and time) in question. ; ; OUTPUTS: ; Ra - The right ascension of the moon at that date in radians. ; Dec - The declination of the moon at that date in radians. ; X, Y, Z - the geocentric equatorial inertial coordinates of the ; center of the Moon w.r.t. the center of the Earth (Re). ; ; SIDE EFFECTS: ; The user should be aware that the equations used are not highly ; accurate. The error in right ascension will rarely exceed 0.3 ; degrees; the error in declination will rarely exceed 0.2 degrees. ; ; PROCEDURE: ; The lunar position is computed from the low precision formulae provided ; in the Astronomical Almanac. ; ; EXAMPLE: ; IDL> juldate,[1982,4,6],jd ;Get reduced julian date ; IDL> jd = jd + 2400000. ;Convert to full julian date ; IDL> moonpos, jd, ra ,dec ;Get RA and Dec of moon ; IDL> print,adstring(ra*!RADEG,dec*!RADEG) ; ==> 11 17 6.1 +09 17 56.0 ; ; This is about 5' from the position given in the Astronomical Almanac ; MODIFICATION HISTORY: ; Written by Michael R. Greason, STX, 31 October 1988. ; Modified by David P. Steele, CNSR, 17 February 1994: ; Added code to calculate and return the GEI cartesian ; coordinates of the Moon, using formulae from the Astronomical ; Almanac. ;- On_error,2 if N_params() EQ 0 then begin print,'Syntax - MOONPOS, dte, ra, dec' print,'dte = Julian Date, Output Ra and Dec are in radians' return endif t = (double(dte) - 2451545.0D0) / 36525.D0 d2r = !DPI / 180.D0 sz = size(dte) scalar = sz(0) eq 0 ;Scalar input? if scalar then t = dblarr(1) + t ; theta1 = (134.9D0 + (477198.85D0 * t)) * d2r lamda = 218.32D0 + (481267.883D0 * t) + (6.29D0 * sin(theta1)) theta1 = (259.2D0 - (413335.38D0 * t)) * d2r theta2 = (235.7D0 + (890534.23D0 * t)) * d2r lamda = lamda - (1.27D0 * sin(theta1)) + (0.66D0 * sin(theta2)) theta1 = (269.9D0 + (954397.70D0 * t)) * d2r theta2 = (357.5D0 + (35999.05D0 * t)) * d2r lamda = lamda + (0.21D0 * sin(theta1)) - (0.19D0 * sin(theta2)) theta1 = (186.6D0 + (966404.05D0 * t)) * d2r lamda = lamda - (0.11D0 * sin(theta1)) neg = where( lamda lt 0.D0, Nneg) if ( Nneg GT 0 ) then lamda(neg) = 360.D0 + ( lamda(neg) mod 360) lamda = (lamda mod 360.D0) * d2r ; theta1 = (93.3D0 + (483202.03D0 * t)) * d2r theta2 = (228.2D0 + (960400.87D0 * t)) * d2r beta = (5.13D0 * sin(theta1)) + (0.28D0 * sin(theta2)) theta1 = (318.3D0 + (6003.18D0 * t)) * d2r theta2 = (217.6D0 - (407332.20D0 * t)) * d2r beta = beta - (0.28D0 * sin(theta1)) - (0.17D0 * sin(theta2)) neg = where( beta LT 0.D0, Nneg ) if ( Nneg GT 0 ) then beta(neg) = 360.D0 + ( beta(neg) mod 360) beta = (beta mod 360.D0) * d2r theta1 = (134.9 + 477198.85D0 * t)*d2r theta2 = (259.2 - 413335.38D0 * t)*d2r pi = (0.9508 + 0.0518*COS(theta1) + 0.0095*COS(theta2)) theta1 = (235.7 + 890534.23D0 * t)*d2r theta2 = (269.9 + 954397.70D0 * t)*d2r pi = pi + 0.0078*COS(theta1) + 0.0028*COS(theta2) r = 1.0D0/SIN(pi*d2r) l = cos(beta) * cos(lamda) m = (0.9175D0 * cos(beta) * sin(lamda)) - (0.3978D0 * sin(beta)) n = (0.3978D0 * cos(beta) * sin(lamda)) + (0.9175D0 * sin(beta)) x = r*l y = r*m z = r*n ra = atan(m, l) neg = where( ra LT 0.0D, Nneg ) if ( Nneg GT 0 ) then ra(neg) = 2.D0*!DPI + ra(neg) dec = asin(n) if scalar then begin ra = ra(0) & dec = dec(0) & x = x(0) & y = y(0) & z = z(0) endif return end ; -------------------------------------------------------------------- ;+ ; NAME: ; LMST ; PURPOSE: ; Give local mean sidereal time. ; CATEGORY: ; CALLING SEQUENCE: ; st = lmst( jd, ut, long) ; INPUTS: ; jd = Julian Day (starting at noon). in ; ut = Universal time as fraction of day. in ; long = observer longitude (deg, East is +). in ; KEYWORD PARAMETERS: ; OUTPUTS: ; st = sidereal time as fraction of day. out ; COMMON BLOCKS: ; NOTES: ; MODIFICATION HISTORY: ; R. Sterner. 16 Sep, 1986. ; R. Sterner, 15 Jan, 1991 --- converted to V2. ; Johns Hopkins University Applied Physics Laboratory. ; ; Copyright (C) 1986, Johns Hopkins University/Applied Physics Laboratory ; This software may be used, copied, or redistributed as long as it is not ; sold and this copyright notice is reproduced on each copy made. This ; routine is provided as is without any express or implied warranties ; whatsoever. Other limitations apply as described in the file disclaimer.txt. ;- ;------------------------------------------------------------- FUNCTION LMST, JD, UT, LNG, help=hlp if (n_params(0) lt 3) or keyword_set(hlp) then begin print,' Give local mean sidereal time.' print,' st = lmst( jd, ut, long)' print,' jd = Julian Day (starting at noon). in' print,' ut = Universal time as fraction of day. in' print,' long = observer longitude (deg, East is +). in' print,' st = sidereal time as fraction of day. out' return, -1 endif SOL2SID = 1.0027379093D0 ; solar to sidereal rate. T0 = 6D0/24. + 38D0/1440. + 45.836D0/86400. ; 0 pt sidereal time. DLONG = -LNG/360.D0 ; time diff from Greenwich. TDAYS = JD - 2415020.5 ; days since noon Jan 0, 1986. T = T0 + (SOL2SID - 1.0)*TDAYS + SOL2SID*UT - DLONG RETURN, T - FLOOR(T) ; only want fraction. END ; -------------------------------------------------------------------- ;+ ; NAME: ; MOONAZEL ; PURPOSE: ; Calculation of the direction to the Moon at a given time and ; place. ; CATEGORY: ; ; CALLING SEQUENCE: ; MOONAZEL, YEAR, MONTH, DAY, HOUR, MINUTE, SECOND $ ; , AZIMUTH, ELEVATION [, MRA, MDEC] ; INPUTS: ; YEAR: year, in either 4-digit or 2-digit representation ; MONTH: month, a number from 1 to 12 ; DAY: day of month ; HOUR: hour, in Universal Time ; MINUTE: minute, in UT ; SECOND: second, in UT ; OPTIONAL INPUTS: ; None. ; OUTPUTS: ; AZIMUTH: The angle East of geographic North of the lunar ; meridian ; ELEVATION: The angle up from the horizon to the Moon in the ; meridian specified by AZIMUTH (ELEVATION < 90 deg) ; OPTIONAL OUTPUTS: ; SRA: The lunar right ascension ; SDEC: The lunar declination ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; None. ; RESTRICTIONS: ; The algorithm giving the lunar position is valid for an ; uncertain period of time. The accuracy of the position is ; typically better than 0.3 deg in RA, and 0.2 deg in ; declination. The place for which the lunar direction is ; calculated is assumed to be the location of the Polar ; Camera, and this is inferred from the specified date. ; PROCEDURE: ; The location of the Polar Camera is established from the ; date specified. The direction to the Moon in inertial ; coordinates is calculated using a routine coded by Michael ; R. Greason, Hughes STX Corporation, October 1988. The ; rest of the calculation is based on equations in the ; Observer's Handbook. ; EXAMPLE: ; moonazel,94,10,12,0,0,0,az,el,/ref ; calculate lunar ; ; az and el for ; ; October 12, 1994, ; ; 00:00 UT (Eureka). ; SEE ALSO: ; moonpos.pro ; MODIFICATION HISTORY: ; Written by: D P Steele, ISR, U of Calgary, March 1994. ; Documented: September 1994. ;- PRO moonazel,year,month,day,hour,minute,second,azimuth,elevation $ ,mra,mdec ; This routine computes the APPROXIMATE lunar azimuth and elevation ; for a given observing location. If enough arguments are specified, ; the lunar right ascension and declination returned by moonpos.pro ; are also returned. fhour=hour+(second/60.+minute)/60. yr=year l1900=WHERE(yr LT 1900.,nl1900) IF nl1900 GT 0 THEN BEGIN l50=WHERE(yr(l1900) LT 50.,nl50) IF nl50 GT 0 THEN yr(l1900(l50))=yr(l1900(l50))+2000. IF nl50 LT nl1900 THEN BEGIN g50=setminus(INDGEN(N_ELEMENTS(l1900)),l50) yr(l1900(g50))=yr(l1900(g50))+1900. ENDIF ENDIF PCsite=PoCasite(year,month,day) ns=N_ELEMENTS(PCsite) lat=REPLICATE(51.08D0,ns) ; default to Calgary lon=REPLICATE(245.8681D0,ns) rabb=WHERE(PCsite EQ 2,nrabb) ; Rabbit Lake coordinates IF nrabb GT 0 THEN BEGIN lat(rabb)=58.23D0 lon(rabb)=256.333D0 ENDIF eur=WHERE(PCsite EQ 3,neur) ; Eureka coordinates IF neur GT 0 THEN BEGIN lat(eur)=80.0533D0 lon(eur)=273.5839D0 ENDIF other=WHERE(ABS(PCsite-2) GT 1,nother) IF nother GT 0 THEN BEGIN PRINT,'An unrecognized Polar Camera site is indicated!' PRINT,'Enter the N lat and E long (deg)' READ,olat,olon lat(other)=olat lon(other)=olon END IF ns EQ 1 THEN BEGIN ; scalar inputs, scalar outputs lat=lat(0) lon=lon(0) ENDIF jdf=ymdh2jdf(yr,month,day,fhour) moonpos,jdf,mra,mdec,mx,my,mz ; mra, mdec in radians; mx, my, mz in Re t2000=(DOUBLE(jdf)-2451545.0D0)/36525.0D0 locmeanst=lmst(ROUND(jdf),fhour/24.0D0,lon)*360.D0 ; lmst in degrees xp=mx-COS(lat*!DTOR)*COS(locmeanst*!DTOR) ; topocentric yp=my-COS(lat*!DTOR)*SIN(locmeanst*!DTOR) ; EI coordinates zp=mz-SIN(lat*!DTOR) ; wrt lat, lon rp=SQRT(xp*xp+yp*yp+zp*zp) alphap=ATAN(yp,xp) ; apparent RA wneg=WHERE(alphap LT 0,nwneg) IF nwneg GT 0 THEN alphap(wneg)=alphap(wneg)+2.0*!DPI deltap=ASIN(zp/rp) ; apparent dec mra=alphap mdec=deltap hourangle=(locmeanst-!RADEG*mra) MOD 360.0D0 ; ditto for hourangle whan=WHERE(hourangle LT 0.0D0,nwhan) IF nwhan GT 0 THEN hourangle(whan)=hourangle(whan)+360.0D0 elevation=!RADEG*ASIN(SIN(mdec)*SIN(!DTOR*lat) $ +COS(!DTOR*hourangle)*COS(mdec)*COS(!DTOR*lat)) azimuth=!RADEG*ASIN(-COS(mdec)*SIN(!DTOR*hourangle)/COS(!DTOR*elevation)) ca=WHERE((((0.0D0 LE azimuth) AND (azimuth LT 90.0D0)) AND $ ((270.0D0 LE hourangle) AND (hourangle LT 360.0D0))) OR $ (((-90.0D0 LE azimuth) AND (azimuth LT 0.0D0)) AND $ ((0.0D0 LE hourangle) AND (hourangle LT 90.0D0))),nca) IF nca GT 0 THEN azimuth(ca)=180.0D0-azimuth(ca) lz=WHERE(azimuth LT 0.0D0,nlz) IF nlz GT 0 THEN azimuth(lz)=azimuth(lz)+360.0D0 RETURN END ; -------------------------------------------------------------------- ;+ ; Name: ; legend ; Purpose: ; This procedure makes a legend for a plot. The legend can contain ; a mixture of symbols, linestyles, and filled polygons (usersym). ; Examples: ; The call: ; legend,['diamond','asterisk','square'],psym=[4,2,6] ; produces: ; ----------------- ; | | ; | <> diamond | ; | * asterisk | ; | [] square | ; | | ; ----------------- ; Each symbol is drawn with a plots command, so they look OK. ; Other examples are given in usage. ; Usage: ; legend,items,linestyle=linestyle ; vertical legend at upper left ; legend,items,psym=psym ; ditto except using symbols ; legend,items,psym=psym,/horizontal ; horizontal format ; legend,items,psym=psym,box=0 ; sans border ; legend,items,psym=psym,delimiter='=' ; embed an '=' betw psym & text ; legend,items,psym=psym,margin=2 ; 2-character margin ; legend,items,psym=psym,position=pos ; position of legend ; legend,items,psym=psym,number=2 ; plot two symbols, not one ; legend,items,/fill,psym=[8,8,8],colors=[10,20,30]; 3 filled squares ; Inputs: ; items = text for the items in the legend, a string array ; Optional Inputs: ; linestyle = array of linestyle numbers ; psym = array of plot symbol numbers. If psym(i) is negative, then a ; line connects pts for ith item. If psym(i) = 8, then the ; procedure usersym is called with vertices define in the ; keyword usersym. ; N. B.: Choose either linestyle, psym or both. If both linestyle and ; psym parameters are present, normal plot behaviour occurs. ; By default, if psym is positive, you get one point so there is ; no connecting line. ; Optional Keywords: ; /help = flag to print header ; /horizontal = flag to make the legend horizontal ; /vertical = flag to make the legend vertical (D=vertical) ; box = flag to include/omit box around the legend (D=include) ; delimiter = embedded character(s) between symbol and text (D=none) ; colors = array of colors for plot symbols/lines (D=!color) ; textcolors = array of colors for text (D=!color) ; margin = margin around text measured in characters and lines ; spacing = line spacing (D=bit more than character height) ; pspacing = psym spacing (D=2 characters) ; charsize = just like !p.charsize for plot labels ; position = normalized coordinates of the upper left of the legend ; number = number of plot symbols to plot or length of line (D=1) ; usersym = 2-D array of vertices, cf. usersym in IDL manual. (D=square) ; /fill = flag to fill the usersym ; Outputs: ; legend to current plot device ; Common blocks: ; none ; Procedure: ; If keyword help is set, call doc_library to print header. ; Restrictions: ; Here are some things that aren't implemented. ; It would be nice to allow data and device coords as well. ; An orientation keyword would allow lines at angles in the legend. ; An array of usersyms would be nice---simple change. ; An order option to interchange symbols and text might be nice. ; Somebody might like double boxes, e.g., with box = 2. ; Side Effects: ; Modification history: ; write, 24-25 Aug 92, F K Knight (knight@ll.mit.edu) ;- pro legend,help=help,items,linestyle=linestyle,psym=psym $ ,horizontal=horizontal,vertical=vertical,box=box,margin=margin $ ,delimiter=delimiter,spacing=spacing,charsize=charsize,pspacing=pspacing $ ,position=position,number=number,colors=colors,textcolors=textcolors $ ,fill=fill,usersym=usersym ; ; =====>> HELP ; if keyword_set(help) then begin & doc_library,'legend' & return & endif ; ; =====>> SET DEFAULTS ; n = n_elements(items) szt = size(items) if (szt(szt(0)+1) ne 7) then message,'First parameter must be a string array. For help, type legend,/help' doline = (n_elements(psym) ne n) if (n_elements(linestyle) ne n) and doline then message,'You must have either linestyle or psym keyword with the same number of elements as items. For help, type legend,/help' if n_elements(linestyle) eq 0 then linestyle = intarr(n) ; D=SOLID if n_elements(psym) eq 0 then psym = intarr(n) ; D=SOLID if n_elements(horizontal) eq 0 then begin ; D=VERTICAL if n_elements(vertical) eq 0 then vertical = 1 endif else begin if n_elements(vertical) eq 0 then vertical = not horizontal endelse if n_elements(box) eq 0 then box = 1 if n_elements(margin) eq 0 then margin = 1 if n_elements(delimiter) eq 0 then delimiter = '' if n_elements(charsize) eq 0 then charsize = 1 if n_elements(spacing) eq 0 then spacing = 1.2 if n_elements(pspacing) eq 0 then pspacing = 3 xspacing = !d.x_ch_size/float(!d.x_size) * (spacing > charsize) yspacing = !d.y_ch_size/float(!d.y_size) * (spacing > charsize) if !x.window(0) eq !x.window(1) then begin plot,/nodata,xstyle=4,ystyle=4,[0],/noerase endif ; next line takes care of weirdness with small windows pos = [min(!x.window),min(!y.window),max(!x.window),max(!y.window)] if n_elements(position) eq 0 then position = [pos(0),pos(3)] + [xspacing,-yspacing] if n_elements(number) eq 0 then number = 1 if n_elements(colors) eq 0 then colors = !color + intarr(n) if n_elements(textcolors) eq 0 then textcolors = !color + intarr(n) fill = keyword_set(fill) if n_elements(usersym) eq 0 then usersym = [[0,0],[0,1],[1,1],[1,0]]-0.5 ; ; =====>> INITIALIZE POSITIONS ; yoff = 0.2*yspacing maxwidth = 0 ; SAVED WIDTH FOR DRAWING BOX x0 = position(0) + (margin)*xspacing ; INITIAL X & Y POSITIONS y0 = position(1) - (margin-0.5)*yspacing y = y0 ; STARTING X & Y POSITIONS x = x0 if vertical then begin ; CALC OFFSET FOR TEXT START xt = 0 ; DEFAULT X VALUE for i = 0,n-1 do begin if psym(i) eq 0 then num = (number + 1) > 3 else num = number if psym(i) lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE if psym(i) eq 0 then expand = 1 else expand = 2 xt = (expand*pspacing*(num-1)*xspacing) > xt endfor endif ; NOW xt IS AN X OFFSET TO ALIGN ALL TEXT ENTRIES ; ; =====>> OUTPUT TEXT FOR LEGEND, ITEM BY ITEM ; =====>> FOR EACH ITEM, PLACE SYM/LINE, THEN DELIMITER, ; =====>> THEN TEXT---UPDATING X & Y POSITIONS EACH TIME. ; for i = 0,n-1 do begin if vertical then x = x0 else y = y0 ; RESET EITHER X OR Y x = x + xspacing y = y - yspacing if psym(i) eq 0 then num = (number + 1) > 3 else num = number if psym(i) lt 0 then num = number > 2 ; TO SHOW CONNECTING LINE if psym(i) eq 0 then expand = 1 else expand = 2 xp = x + expand*pspacing*indgen(num)*xspacing yp = y + yoff + intarr(num) if psym(i) eq 0 then begin xp = [min(xp),max(xp)] ; TO EXPOSE LINESTYLES yp = [min(yp),max(yp)] ; DITTO endif if psym(i) eq 8 then usersym,usersym*charsize,fill=fill,color=colors(i) plots,xp,yp,color=colors(i),/normal $ ,linestyle=linestyle(i),psym=psym(i),symsize=charsize if vertical then x = x + xt else x = max(xp) x = x + xspacing xyouts,x,y,delimiter,width=width,/norm,color=textcolors(i),size=charsize x = x + width + xspacing xyouts,x,y,items(i),width=width,/norm,color=textcolors(i),size=charsize x = x + width if not vertical and (i lt (n-1)) then x = x+2*xspacing; ADD INTER-ITEM SPACE maxwidth = (x + xspacing*margin) > maxwidth ; UPDATE MAXIMUM X endfor ; ; =====>> OUTPUT BORDER ; if box then begin x = position(0) y = position(1) if vertical then bottom = n else bottom = 1 ywidth = - (2*margin+bottom-0.5)*yspacing plots,[x,maxwidth,maxwidth,x,x],y + [0,0,ywidth,ywidth,0],/norm endif return end ; -------------------------------------------------------------------- ;+ ; NAME: ; PSOPEN ; PURPOSE: ; Opens a file for PostScript graphics output ; CATEGORY: ; Utility ; CALLING SEQUENCE: ; PSOPEN,pfile,font=font,portrait=portrait,landscape=landscape $ ; ,color=color,longpage=longpage ; INPUTS: ; pfile = string specifying name of PostScript output file to create ; KEYWORD PARAMETERS: ; font: specifies PostScript font other than default (e.g., 'times') ; portrait: specifies portrait orientation on page (default) ; landscape: specifies landscape orientation on page ; color: specifies color PostScript output, 8 bits per pixel ; longpage: specifies full length of page available for plotting ; ; Available fonts include AVANTGARDE, BKMAN, COURIER, HELVETICA, ; PALATINO, SCHOOLBOOK, TIMES, ZAPFCHANCERY, ZAPFDINGBATS ; ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; COMMON BLOCKS: ; PS_STUFF,Pfont,Dname,Ps_filename ; used by CLOSE_PS to reset font, graphics device type, and to enable ; it to print a ^D to the file for compatibility with QMS printer ; SIDE EFFECTS: ; Opens a file for plotting; sets plot type to PostScript; ; sets !P.FONT to 0 (if a font change requested) ; RESTRICTIONS: ; None known. ; PROCEDURE: ; ; MODIFICATION HISTORY: ; T.J.Hughes, HIA/NRCC , November 30, 1990 ; D. P. Steele, CNSR/UofC, April 21 1993, modified to accept command ; line parameters; added color output option ;- PRO PSOPEN,pfile,font=font,portrait=portrait,landscape=landscape $ ,color=color,longpage=longpage COMMON ps_stuff,pfont,dname,ps_filename nfont='' IF KEYWORD_SET(font) THEN nfont=',/'+font orient=',/PORTRAIT' IF KEYWORD_SET(landscape) THEN orient=',/LANDSCAPE' ystuff='' IF KEYWORD_SET(longpage) THEN ystuff=',YSIZE=22.6,YOFFSET=3.46' ps_filename=pfile dvfile="'"+pfile+"'" dname=!D.NAME pfont=!P.FONT SET_PLOT,'ps' dev_string='DEVICE,FILE='+dvfile+orient+nfont+ystuff IF KEYWORD_SET(color) THEN dev_string=dev_string+',/color,bits_per_pixel=8' istat=EXECUTE(dev_string) RETURN END ; -------------------------------------------------------------------- PRO PSCLOSE, NORMPS3=NORMPS3 ;+ ; NAME: ; PSCLOSE ; PURPOSE: ; Closes a PostScript output file (opened by OPEN_PS). It then reopens ; the file for Update, and appends a D to the file for compatability ; with the QMS PostSCript printer ; CATEGORY: ; Utility ; CALLING SEQUENCE: ; PSCLOSE ; INPUTS: ; None ; OPTIONAL INPUT PARAMETERS: ; None ; KEYWORD PARAMETERS: ; /NORMPS3: If set, prevents the output PostScript file being ; reopened to remove the spurious PS-Adobe-3.0 ; characters in the first line that mess up the ; PostScript printers. ; OUTPUTS: ; None ; OPTIONAL OUTPUT PARAMETERS: ; None ; COMMON BLOCKS: ; PS_STUFF,pfont,dname,ps_filename ; SIDE EFFECTS: ; Closes the file PS_FILENAME ; Appends a D to that file ; RESTRICTIONS: ; None known ; PROCEDURE: ; ; MODIFICATION HISTORY: ; T.J.Hughes, HIA/NRCC , November 30, 1990 ; D.P.Steele, CNSR/UofC, April 21, 1993, modified from CLOSE_PS.PRO ;- COMMON ps_stuff,pfont,dname,ps_filename ; +-------------------------------------------------------------------+ ; |Close PS_FILENAME so that all the graphics buffers are cleared. | ; |EMPTY didn't seem to do it, and we can't use FLUSH,UNIT because | ; |we don't know what unit has been used. | ; +-------------------------------------------------------------------+ DEVICE,/CLOSE ctrl_d=string(4b) ; +-------------------------------------------------------------------+ ; |Open the same file again, for updating, and append a D to | ; |it | ; +-------------------------------------------------------------------+ GET_LUN,unit OPENU,unit,ps_filename,/APPEND PRINTF,unit,ctrl_d CLOSE,unit FREE_LUN,unit ; +-------------------------------------------------------------------+ ; |Reset the output graphics type to what it was before OPEN_PS was | ; |called, and reset !P.FONT similarly | ; +-------------------------------------------------------------------+ SET_PLOT,dname !P.FONT=pfont ; +-------------------------------------------------------------------+ ; |Edit the first line of the file to remove the PS-Adobe-3.0 string | ; |that messes up (at least) the b/w PostScript printer. | ; +-------------------------------------------------------------------+ IF KEYWORD_SET(normps3) THEN GOTO,done pstext=qgetfile(ps_filename) pstext(0)=STRMID(pstext(0),0,2) mputfile,ps_filename,pstext done: RETURN END ; -------------------------------------------------------------------- PRO qlook,file,debug=debug,ps=ps ; IDL Version 3.1.0 (RISC/os mipseb) ; Journal File for cnsr3@jasper ; Working directory: /jasper/cnsr3_data1 ; Date: Mon Jan 9 10:18:08 1995 ; First convert HH:MM:SS times to UT seconds. ; Then separate data by camera. ; Then determine cycles and identify dark counts. ; Save dark counts and remove them from the other data. ; Subtract interpolated dark counts from all data. ; Plot data by camera/filter. ; Print usage message IF N_PARAMS() LT 1 THEN BEGIN MESSAGE,"Usage: qlook,file,debug=debug,ps=ps" RETURN ENDIF ; Load the quick-look data file f=rfindfile(file,count=nf) IF nf EQ 0 THEN BEGIN MESSAGE,file+' not found' RETURN ENDIF q=qgetfile(file) nq=N_ELEMENTS(q) IF KEYWORD_SET(debug) THEN help,q day=FIX(STRMID(file,STRLEN(file)-2,2)) READ,'Enter the year (2 or 4 digits) and month number: ',year,month READ,'Enter the CCD pixel binning factor for image pixels (1/2/4): ',binfac ; Extract numeric data and convert HH:MM:SS times to UT seconds cam=INTARR(nq) filt=cam expos=cam data=LONARR(nq) utsec=FLTARR(nq) FOR i=0,nq-1 DO BEGIN qi=q(i) cam(i)=FIX(STRMID(qi,0,1)) filt(i)=FIX(STRMID(qi,2,1)) expos(i)=FIX(STRMID(qi,5,3)) hut=FIX(STRMID(qi,9,2)) mut=FIX(STRMID(qi,12,2)) sut=FIX(STRMID(qi,15,2)) utsec(i)=(hut*60.+mut)*60.+sut data(i)=LONG(STRMID(qi,18,STRLEN(qi)-18)) ENDFOR ; Separate data by camera. c0=WHERE((cam EQ 0) AND (expos NE 0)) nc0=N_ELEMENTS(c0) c1=WHERE((cam EQ 1) AND (expos NE 0)) nc1=N_ELEMENTS(c1) noerr=(nc0 EQ nc1) IF noerr THEN timediff=ABS(utsec(c0)-utsec(c1)) ELSE timediff=2. noerr=noerr AND (MAX(timediff) LE 1) IF NOT noerr THEN BEGIN MESSAGE,'At least one missing image; synch problems' $ ,/INFORMATIONAL STOP RETURN ENDIF ELSE BEGIN npairs=nc0 utsec=utsec(c0) exp0=expos(c0) exp1=expos(c1) f0=filt(c0) f1=filt(c1) data0=data(c0) data1=data(c1) ENDELSE ; Determine cycles and identify dark counts dutsec=(utsec-SHIFT(utsec,1))>0 cycstart=[0,WHERE(dutsec GT 68),npairs] ncycle=N_ELEMENTS(cycstart) nincycle=(SHIFT(cycstart,-1)-cycstart)>0 nincycle=nincycle(0:ncycle-2) cychist=HISTOGRAM(nincycle,MIN=0,BIN=1) nincnorm=MAX(cychist,normlen) longer=WHERE(nincycle GT normlen) dark=cycstart(longer+1)-1 ; indices of dark counts dkut=utsec(dark) ; times of dark counts ; Save dark counts and remove them from the other data dark0=data0(dark) ; values of dark counts dark1=data1(dark) allcts=REPLICATE(1,npairs) allcts(dark)=0 ; 1 if shutter open whimg=WHERE(allcts) utsec=utsec(whimg) exp0=exp0(whimg) exp1=exp1(whimg) f0=f0(whimg) f1=f1(whimg) data0=data0(whimg) ; delete dark counts data1=data1(whimg) ; from data arrays ; Smooth dark counts to remove hot-pixel noise dark0=MEDIAN(dark0,3) dark1=MEDIAN(dark1,3) ; Subtract interpolated dark counts from all data idark0=FIX(interpol(dark0,dkut,utsec)) idark1=FIX(interpol(dark1,dkut,utsec)) data0=(data0-idark0)>0 data1=(data1-idark1)>0 IF KEYWORD_SET(debug) THEN HLP,data0,idark0,data1,idark1 ; Plot data by camera/filter ; Calibration values (Rayleigh-seconds per DN for 2x2 binning) cal0=[129.4, 62.3, 72.4, 62.7, 0.0]*(2./binfac)^2 cal1=[113.3,103.8, 88.0, 69.2, 49.1]*(2./binfac)^2 ; Wavelengths and features fname0=['!17829 nm OH','!17608 nm bg','!17630.0 nm OI' $ ,'!17668 nm N!i2!n','!17520.0 nm NI'] fname1=['!17835 nm OH','!17820 nm bg','!17785 nm N!i2!e+!n' $ ,'!17557.7 nm OI','!17714 nm bg'] ; Prepare plot device IF KEYWORD_SET(ps) THEN BEGIN dot=STRPOS(file,'.') yymm=STRING(year MOD 100,month,FORMAT='("_",2I2.2)') psfile=STRMID(file,0,dot)+yymm+STRMID(file,dot+1,2)+'.ps' psopen,psfile,/long,font='times' ENDIF ELSE WINDOW,0,XSIZE=640,YSIZE=900 ; Set up plot space !P.MULTI=[0,1,2,0,0] date=ymd2date(year,month,day,format='d$ n$ y$') !X.TITLE='!17Universal Time (hours, '+date+')' !X.RANGE=[0,24] !X.STYLE=1 !X.TICKS=4 !X.TICKV=INDGEN(5)*6 !X.MINOR=6 !Y.TITLE='!17Zenith Intensity (R or R/nm)' !Y.RANGE=[0,1000] !Y.STYLE=1 IF KEYWORD_SET(debug) THEN BEGIN STOP resp='' READ,'Enter an alternate full-scale brightness? ',resp IF resp NE '' THEN !Y.RANGE(1)=FLOAT(resp) ENDIF ; Plot data for Camera 0 PLOT,[0,1],/NODATA,TITLE='!17Quick Look Plot for Camera 0, '+date found=REPLICATE(0,5) FOR jf=0,4 DO BEGIN wjf=WHERE(f0 EQ jf,nwjf) IF nwjf GT 0 THEN BEGIN found(jf)=nwjf tf=utsec(wjf) df=data0(wjf)*cal0(jf)/exp0(wjf) IF jf EQ 1 THEN df=df/4.4 ; R => R/nm for 608 nm bg channel IF KEYWORD_SET(debug) THEN STOP OPLOT,tf/3600.,df,LINESTYLE=jf ENDIF ENDFOR ; If any data were plotted, create a legend. IF TOTAL(found) GT 0 THEN BEGIN wfound=WHERE(found) labels=fname0(wfound) lines=wfound syms=REPLICATE(0,N_ELEMENTS(wfound)) ; If Moon up, show its elevation y=REPLICATE(year,25) mo=REPLICATE(month,25) d=REPLICATE(day,25) h=FINDGEN(25) mi=FLTARR(25) s=mi moonazel,y,mo,d,h,mi,s,maz,mel moonup=(MAX(mel) GT 0.) IF moonup THEN BEGIN OPLOT,h,mel*10.,PSYM=-6 labels=[labels,'!17Moon el.*10'] lines=[lines,0] syms=[syms,-6] ENDIF ; Determine the normal coordinates of the upper left plot corner ; and plot the legend xdvpos=!P.CLIP(0) ydvpos=!P.CLIP(3) npos=CONVERT_COORD([xdvpos,ydvpos],/DEVICE,/TO_NORMAL) npos=npos(0:1) legend,labels,linestyle=lines,psym=syms,position=npos ENDIF ; Plot data for Camera 1 PLOT,[0,1],/NODATA,TITLE='!17Quick Look Plot for Camera 1, '+date found=REPLICATE(0,5) FOR jf=0,4 DO BEGIN wjf=WHERE(f1 EQ jf,nwjf) IF nwjf GT 0 THEN BEGIN found(jf)=nwjf tf=utsec(wjf) df=data1(wjf)*cal1(jf)/exp1(wjf) IF jf EQ 1 THEN df=df/8.2 ; R => R/nm for 820 nm bg channel IF jf EQ 4 THEN df=df/6.9 ; R => R/nm for 714 nm bg channel IF KEYWORD_SET(debug) THEN STOP OPLOT,tf/3600.,df,LINESTYLE=jf ENDIF ENDFOR ; If any data were plotted, create a legend. IF TOTAL(found) GT 0 THEN BEGIN wfound=WHERE(found) labels=fname1(wfound) lines=wfound syms=REPLICATE(0,N_ELEMENTS(wfound)) ; Determine the normal coordinates of the upper left plot corner ; and plot the legend xdvpos=!P.CLIP(0) ydvpos=!P.CLIP(3) npos=CONVERT_COORD([xdvpos,ydvpos],/DEVICE,/TO_NORMAL) npos=npos(0:1) legend,labels,linestyle=lines,psym=syms,position=npos ENDIF ; Close the PS device, if necessary. IF KEYWORD_SET(ps) THEN psclose ; We're done RETURN END