;+ ; NAME: ; smei_sky_track ; PURPOSE: ; Track location on sky ; CATEGORY: ; camera/idl/sky ; CALLING SEQUENCE: PRO smei_sky_track, wanted_map , $ body = body , $ degrees = degrees , $ rmbkgnd = rmbkgnd , $ rmzld = rmzld , $ range = range , $ destination = destination , $ equator = equator , $ fix_sun = fix_sun , $ step = step , $ rm_min_count= rm_min_count , $ _extra = _extra ; INPUTS: ; wanted_map scalar, array[1]; type: numerical, time structure, ; or string ; indicates the orbit(s) to be processed ; (passed to smei_getfile together with _extra). ; OPTIONAL INPUT PARAMETERS: ; body=body array[2], type: float ; ecliptic or equatorial coordinates ; of fixed point on sky ; scalar; type: string ; name of one of the bodies for which ; an ephemeris is available from href=big_eph= ; /degrees if set, all angles are in degrees ; (default: radians) ; /rmbkgnd if set, subtract sidereal background ; (for 'sky' maps this is SET by default; ; for 'equ' map this is NOT set by default) ; /rmzld if set, subtract zodiacal light distribution ; by default this is SET for both 'sky' ; and 'equ' maps ; range=range array[2]; type: float; ; default: [-30,390] for 'sky' maps; ; [-10, 10] for 'equ' maps ; brightness range in ADUS ; destination=destination ; name of existing directory. If set this is ; where gif images will be written ; /equator keep horizontal parallel to equatorial plane ; (default: parallel to ecliptic plane) ; /fix_sun if set the direction to the Sun is kept on ; the x-axis (left in a fish-eye plot) ; step=step scalar; type: integer; default: 1 ; step size for processing the selected ; map: only one of every 'step' maps is processed ; rm_min_count=rm_min_count ; scalar; type: integer: default: 0 ; the first 'rm_min_count' maps are used to create ; a 'minimum skymap'. The minimum map is then ; subtracted from the displayed map for all subsequent ; skymaps. ; _extra=_extra ; OUTPUTS: ; OPTIONAL OUTPUT PARAMETERS: ; INCLUDE: @compile_opt.pro ; On error, return to caller ; CALLS: ; InitVar, IsType, ToRadians, ToDegrees, AngleUnits ; smei_sky, smei_getfile, jpl_body, big_eph ; TimeGet, TimeSet, TimeUnit, smei_filename, CheckDir ; SuperArray, smei_orbit_set, smei_orbit_get ; PROCEDURE: ; mei_sky_track,dest=getenv('PWD'),mode='equ',maxelo=10,/deg, $ ; chars=1.2,range=[-8,29.9],xysize=[400,400],nbin=1, /use_flat ; MODIFICATION HISTORY: ; JUL-2008, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ;- InitVar, degrees, /key rpm = ToRadians(degrees=degrees) dpm = ToDegrees(degrees=degrees) InitVar, equator, /key InitVar, fix_sun, /key InitVar, step , 1 InitVar, body, 'C/2004 Q2 (Machholz)' InitVar, rm_min_count, 0 IF IsType(body,/string) THEN BEGIN CASE body OF 'C/2004 Q2 (Machholz)' : BEGIN InitVar, wanted_map, TimeSet(smei=[9490,11930]) END 'V598 Puppis' : BEGIN InitVar, wanted_map, TimeSet('2007_'+['131','265']) body = smei_star_info(name=body, degrees=degrees, $ catalogue='stars_not_subtracted') END ENDCASE ENDIF given_map = smei_getfile(wanted_map, $ count = count , $ _extra = _extra ) IF count EQ 0 THEN $ RETURN save_img = 0 IF IsType(destination,/defined) THEN BEGIN IF NOT CheckDir(destination,/silent) THEN $ RETURN save_img = 1 ENDIF usec = TimeUnit(/sec) orbit = smei_filename(given_map, mode=mode, ccd=ccd) IF ccd[0] THEN BEGIN message, /info, given_map[0] message, 'you seem to be processing a CCD frame !!!' ENDIF orbit = TimeGet(orbit,/smei) orbit = smei_orbit_get(orbit,/round) orbit = orbit[sort(orbit)] orbit = orbit[uniq(orbit)] count = n_elements(orbit) time = TimeSet(smei=smei_orbit_set(orbit,0.5d0)) mode = mode[0] ; Pick up the ephemeris ; If 'body' is numeric it should be a pair of ra,dec ephemeris = IsType(body,/string) CASE ephemeris OF 0: BEGIN ; J2000 (ra,dec) pair pp = body pp = SuperArray(pp,count,/trail) ; [2,count] END 1: BEGIN ; J2000 ephemeris (equatorial) pp = big_eph(time , $ body = body , $ center = jpl_body(/earth,/string), $ degrees = degrees , $ to_equatorial= equator , $ /precess , $ /to_sphere , $ /onebody , $ /silent ) rr = reform(pp[2,*]) ; Distance pp = pp[0:1,*] ; RA/dec [2,count] END ENDCASE ; Get location of sun at time of maps sun = (big_eph(time , $ body = jpl_body(/sun ,/string), $ center = jpl_body(/earth,/string), $ degrees = degrees , $ to_equatorial = equator , $ /precess , $ /to_sphere , $ /onebody , $ /silent ))[0:1,*] ; RA/dec only [2,count] sun *= rpm ; Convert to radians pp *= rpm dphi = sun[0,*]-pp[0,*] gam = atan( sin(dphi)*cos(sun[1,*]) , sin(sun[1,*])*cos(pp[1,*])-cos(sun[1,*])*sin(pp[1,*])*cos(dphi) ) ; Euler angles to put direction to Sun on the positive x-axis ; (left in the fish-eye plot) dabg = [pp[0,*],-pp[1,*],!dpi/2-gam] ; This keeps the x-axis parallel to the equator/ecliptic. IF NOT fix_sun THEN dabg[2,*] = 0 dabg /= rpm ; Back to mystery units sun /= rpm pp /= rpm user_position = [0.99,0.05] user_align = 1.0 zero_point = 0.0/dpm InitVar, rmbkgnd, mode EQ 'sky' InitVar, rmzld , 1 InitVar, range, [-30,390]*(mode EQ 'sky') + [-10,10]*(mode EQ 'equ') min_count = -1 FOR i=0L,count-1,step DO BEGIN found_body = 0 maps = smei_getfile(orbit[i], $ /exact , $ count = nmap , $ _extra = _extra ) ; Extract time at RA/dec pp[*,i] ; Note that pp probably is an ephemeris location calculate ; at the start time of the orbit. FOR j=nmap-1,0,-1 DO BEGIN exten_no = smei_sky_field(maps[j],/orbsecs) time_map = readfits(maps[j],hdr,exten_no=exten_no,/silent) qq = cvsmei( $ from_equatorial = pp[*,i] , $ degrees = degrees , $ /to_map , $ to_mode=smei_sky_field(maps[j],exten_no=exten_no,/cvsmei_mode), $ /silent) qq = round(qq) qq = time_map[qq[0],qq[1]] ; Seconds since TORIGIN ; A valid time indicates that the body is seen is this map found_body = qq NE 0 IF found_body THEN BEGIN qq = TimeOp(/add, $ TimeSet(fxpar(hdr,'TORIGIN')) , $ TimeSet(/diff,qq,usec)) message, /info, 'time update: '+strjoin(TimeGet([time[i],qq],/_ydoy,upto=usec,/roundt),' -> ') time[i] = qq break ENDIF ENDFOR IF found_body THEN BEGIN ; Update ephemeris location for new time IF ephemeris THEN BEGIN qq = big_eph(time[i] , $ body = body , $ center = jpl_body(/earth,/string), $ degrees = degrees , $ to_equatorial= equator , $ /precess , $ /to_sphere , $ /onebody , $ /silent ) rr[ i] = qq[ 2] ; Distance pp[*,i] = qq[0:1] ; RA/dec ENDIF ; Update location of Sun for new time sun[*,i] = (big_eph(time[i] , $ body = jpl_body(/sun ,/string), $ center = jpl_body(/earth,/string), $ degrees = degrees , $ to_equatorial = equator , $ /precess , $ /to_sphere , $ /onebody , $ /silent ))[0:1] ; RA/dec only ; Update the Euler angles for new time sun [*,i] *= rpm ; Convert to radians pp [*,i] *= rpm dabg[*,i] *= rpm dphi = sun[0,i]-pp[0,i] gam = atan( sin(dphi)*cos(sun[1,i]) , sin(sun[1,i])*cos(pp[1,i])-cos(sun[1,i])*sin(pp[1,i])*cos(dphi) ) ; Euler angles to put direction to Sun on the positive x-axis ; (left in the fish-eye plot) dabg[*,i] = [pp[0,i],-pp[1,i],!dpi/2-gam] ; This keeps the x-axis parallel to the equator/ecliptic. IF NOT fix_sun THEN dabg[2,i] = 0 dabg[*,i] /= rpm ; Back to mystery units sun [*,i] /= rpm pp [*,i] /= rpm ENDIF IF save_img THEN $ imgfile = filepath(root=destination, $ 'c0gif_'+TimeGet(time[i],/_ydoy,upto=usec,/roundt)+'.gif') user_string = AngleUnits(from_angle=pp[*,i] , $ degrees = degrees , $ to_sexages=1-equator, $ to_almanac= equator, $ /strings , $ /xyout , $ upto=usec) IF ephemeris THEN $ user_string = [body, user_string, string(rr[i], format='(F7.3)')+' AU'] user_string = strjoin(user_string,' ') mkplot = found_body AND (rm_min_count EQ 0 OR min_count GE rm_min_count) smei_sky, orbit[i] , $ noplot = 1-mkplot , $ /exact , $ /fixgain , $ /naked , $ /legend , $ equator = equator , $ rmzld = rmzld , $ rmbkgnd = rmbkgnd , $ degrees = degrees , $ dabg = dabg[*,i] , $ zero_point = zero_point , $ imgfile = imgfile , $ range = range , $ time = time[i] , $ user_position = user_position , $ user_align = user_align , $ user_string = user_string , $ skymap = skymap , $ min_map = min_map , $ _extra = _extra IF rm_min_count GT 0 THEN BEGIN IF IsType(min_sav,/undefined) THEN $ min_sav = make_array( $ type = IsType(skymap) , $ dim = [size(skymap,/dim),rm_min_count], $ value = BadValue(skymap)) min_count += 1 min_sav[*,*,min_count mod rm_min_count] = skymap min_map = min(min_sav,3,/nan) ENDIF ENDFOR RETURN & END