PRO yopos, timein, latitude, longitude, p0long=p0long, replot=replot, $ giffile=giffile, zbuffer=zbuffer, outimg=outimg, qdebug=qdebug, $ plotnight=plotnight ; Coefficients for Equation of Time ceot = [-0.120553, -0.0323706, 0.000394917, 0.163718] ; Keplers constant for Earth ; (conversion between orbit period and orbit semi-major axis) kepler = 0.01877 ; Earth mean radius (km) earth_radius = 6371.0 ; Station Code for Kagoshima stcode = BYTE(85) ; Latitude, longitude, of Kagoshima (degrees) ; Latitude = orbital inclination kaglat = 31.345 kaglong = 131. ; Name of world map file (file location is version dependent) worldname = 'worldelv.dat' subdir=(['images',concat_dir('examples','data')])(since_version('4')) worldfile=(concat_dir(!dir,concat_dir(subdir,worldname)))(0) if not file_exist(worldfile) then begin message,/info,"Cannot find elevation file: " + worldname return endif ; Common to keep track of display window COMMON YOPOS_WINDOW, winnumber, startx, starty, newworld ; Scale factors for Yohkoh spacecraft image yoscale = 0.5 ylatscale = 512./(512. - 53.*yoscale) ylongscale = 512./(512. - 97.*yoscale) ; HANDLE TIME ----------------------------------------------------------- ; Is there an input time? IF( N_ELEMENTS(timein) GT 0 ) THEN BEGIN timstr = ANYTIM2INTS( timein ) IF( KEYWORD_SET(replot) EQ 0 ) THEN timestr = timein ENDIF ELSE BEGIN timstr = ANYTIM2INTS( UT_TIME() ) ENDELSE ; Look +/- 12 hours ; slf - fix problem (day boundry) ; GET YOHKOH ORBIT PARAMETERS --------------------------------------------- ; Inclination never changes inclin = kaglat RD_WEEK_FILE, timegrid(timstr,hour=-12),timegrid(timstr,hour=12),'FEM', fem IF( GET_NBYTES(fem) LT 10) THEN BEGIN PRINT,' No FEM data available for this time.' RETURN END nfem = N_ELEMENTS(fem) - 1 ; Sidereal period in minutes period = ( 1440. * (fem(nfem).day - fem(0).day) + $ (fem(nfem).time - fem(0).time)/60000. ) / $ (fem(nfem).sc_rev - fem(0).sc_rev) ; Synodic period persynod = 1./(1./period - 1./1440.) perrat = period/persynod ; Orbit semimajor axis semimajor = SQRT( period^3/kepler ) ; When was Yohkoh over Kagoshima (U.T. seconds)? nstation = N_ELEMENTS(fem(0).st_station) timsec = fem.time/1000. + (fem.day - fem(0).day) * 86400. timsec = TRANSPOSE( REBIN( timsec, nfem+1, nstation) ) contact = WHERE( fem.st$station EQ stcode ) starts = fem.st_station + timsec ends = fem.en_station + timsec duration = ends(contact) - starts(contact) centers = (starts(contact) + ends(contact))/2. ldur = MAX(duration, top) timref = centers(top) ; Choose map central longitude IF (n_elements(p0long) GT 0) then $ plong = p0long $ else $ plong = 139. ; ISAS longitude ; TOP OF REPEAT LOOP REPEATLOOP: ; Force the size scale yos = yoscale ; If time is input, translate, otherwise, use the clock IF (N_ELEMENTS(timestr) EQ 0) THEN BEGIN timenow = UT_TIME() ENDIF ELSE BEGIN timenow = timestr ENDELSE timstr = ANYTIM2INTS( timenow ) uthours = timstr.time/3600000. ptitle = FMT_TIM( timenow ) + ' U.T.' ; MAKE WINDOW ------------------------------------------------------------ ; Make a 512 x 512 window so we can use YO_LOGO later ; Reuse window if it exists and is the right size ; optionally use Z-buffer (via wdef call) zbuffer = keyword_set(zbuffer) or keyword_set(giffile) IF( N_ELEMENTS(winnumber) GT 0 and (1-zbuffer)) THEN BEGIN DEVICE, WINDOW=win_arr IF(win_arr(winnumber)) THEN BEGIN ; "old" window exists wtemp=!d.window wset,winnumber if (!d.x_size ne 512) or (!d.y_size ne 512) then wset,wtemp ENDIF if !d.window ne winnumber then delvarx,winnumber ENDIF if zbuffer then set_plot,'z' if n_elements(winnumber) eq 0 or zbuffer then WDEF, winnumber, 512, 512, zbuffer=zbuffer ; PLOT MAP ------------------------------------------------------------ ; Set the scaling left = plong - 180. right = plong + 180. LOADCT, 8 TEK_COLOR MAP_SET, 0., plong, /CYLINDRICAL, TITLE = ptitle, $ COLOR=1 ; Get the background image IF( N_ELEMENTS(newworld) LE 0) THEN BEGIN world = BYTARR(360,360) OPENR, lun, worldfile, /GET_LUN READU, lun, world CLOSE, lun FREE_LUN, lun sea = WHERE(world LT 125) world(sea) = 4 south = world(*,0:64) hi = WHERE(south ge 125) south(hi) = 255 world(*,0:64) = south north = world(250:*,295:*) hi = WHERE(north GE 140) north(hi) = 255 world(250:*,295:*) = north ; Remap image if since_version('4.') then begin newworld = map_patch(world, $ xstart=startx,ystart=starty, $ lon0=0., lon1=360.) message, /info, 'Inside since_version loop' endif else $ newworld = MAP_IMAGE(world, startx, starty,COMPRESS=1,LONMIN=0., $ /WHOLE_MAP) ENDIF world = 0 TV, newworld, startx, starty MAP_CONTINENTS, COLOR=9 ; PLOT TERMINATOR, DAY/NIGHT ------------------------------------------- ; Call for ephemeris ephsun = GET_SUN( timstr ) rasun = ephsun(6) ras = rasun * !PI / 12. ; Equation of time (in hours) eot = ceot(0)*COS(ras) + ceot(1)*SIN(ras) + ceot(2)*COS(2.*ras) $ + ceot(3)*SIN(2.*ras) ; Longitude of noon zenlong = (12. - (uthours + eot) ) * 15. ; Declination of Sun decsun = ephsun(9) ; Terminator ; Initial coordinate system has Sun in +z direction, ; terminator lies in (x,y) plane. ; Longitude of noon defines zero angle in (x,y) plane ; Rotate by (90 - decsun) around the y axis theta = FINDGEN(360.) * !DTOR x = COS(theta) y = SIN(theta) phi = (90. - decsun) * !DTOR xp = x * COS(phi) yp = y zp = -x * SIN(phi) ; Convert to (latitude, longitude) latitude = ASIN(zp) / !DTOR longitude = ATAN(yp,xp) / !DTOR longitude = longitude + zenlong ; Eliminate discontinuity at theta = 180. ; (Changing the (+) to a (-) fixed some day/night problems - BNH) longitude(180:*) = longitude(180:*) + 360. ; Find plot range, fit longitudes into it ; Since longitudes are in order, just do circular shift over = WHERE(longitude GT right, ocount) IF ( ocount GT 0 ) THEN BEGIN longitude(over) = longitude(over) - 360. longitude = SHIFT(longitude, ocount) latitude = SHIFT(latitude, ocount) ENDIF under = WHERE(longitude LT left, ocount) IF(ocount GT 0) THEN longitude(under) = longitude(under) + 360. ; Pad out the array to get complete plot longt = [longitude, longitude(0)] lat = [latitude, latitude(0)] ; Plot terminator OPLOT, longt, lat, COLOR=0, THICK=2 ; Mark which side is night subsl = SORT(longitude) nlats = N_ELEMENTS(longitude) IF(decsun GE 0.) THEN edge = -90. ELSE edge = 90. if keyword_set(qdebug) then stop ; ------ This Just Doesn't Work. -- (BNH) if keyword_set(plotnight) then $ polyfill, [left,left,longitude(subsl),right,right], $ [edge,latitude(subsl(0)),latitude(subsl),latitude(subsl(nlats-1)),edge], $ /LINE_FILL, COLOR=0, THICK=1, /DATA if keyword_set(qdebug) then stop ; PLOT SOUTH ATLANTIC ANOMALY ----------------------------------------- ; Position, extent are: xsaapos = -25. ysaapos = -30. xsaaxtent = 75. ysaaxtent = 25. theta = FINDGEN(25) * 15. * !DTOR xsaa = xsaapos + xsaaxtent * COS(theta) ysaa = ysaapos + ysaaxtent * SIN(theta) ; Handle spanning the map edge higher = WHERE(xsaa GE left) lower = WHERE(xsaa LE left) IF( (lower(0) LT 0) OR (higher(0) LT 0) ) THEN BEGIN POLYFILL, xsaa, ysaa, /LINE_FILL, COLOR=2, /DATA ENDIF ELSE BEGIN POLYFILL, xsaa(higher), ysaa(higher), /LINE_FILL, COLOR=2, /DATA POLYFILL, [left,left,MIN(xsaa(higher)),MIN(xsaa(higher))], $ [MAX(ysaa(higher)),MIN(ysaa(higher)),MIN(ysaa(higher)), $ MAX(ysaa(higher))], /LINE_FILL, COLOR=2, /DATA POLYFILL, [xsaa(lower),right,right], [ysaa(lower),MIN(ysaa(lower)), $ MAX(ysaa(lower))], /LINE_FILL, COLOR=2, /DATA ENDELSE ; LONGITUDE, LATITUDE LINES MAP_GRID ; PLOT RADIUS OF KSC CONTACTS ----------------------------------------------- ; Use 10 degree elevation angle anga = !DTOR * (90. + 10.) angc = ASIN( earth_radius * SIN(anga)/semimajor ) angb = (!PI - anga - angc)/!DTOR x = kaglong + angb * COS(theta)/COS(kaglat*!DTOR) y = kaglat + angb * SIN(theta) OPLOT, x, y, COLOR=7 ; PLOT YOHKOH ORBIT AND LOCATION ---------------------------------------- ; Time difference in minutes between map time and S/C overhead at KSC diftim = (timstr.time/1000. + (timstr.day - fem(0).day)*86400. - timref)/60. ; Orbital longitude beyond integral number of revolutions revs = 2. * !PI * ( (diftim / period) MOD 1. ) ; Convert orbital longitude to Earth xr = SIN(revs) yr = COS(revs) * COS(inclin) rr = SQRT(xr^2 + yr^2) erevs = ACOS(yr/rr)/!DTOR IF( xr LT 0.) THEN erevs = 360. - erevs ; Earth longitude rotation rotearth = 360. * diftim/1440. ; Spacecraft Earth longitude sclong = (kaglong + erevs - rotearth) MOD 360. IF(sclong LT left) THEN sclong = sclong + 360. IF(sclong GT right) THEN sclong = sclong - 360. orblong = (kaglong - rotearth) MOD 360. IF(orblong LT left) THEN orblong = orblong + 360. IF(orblong GT right) THEN orblong = orblong - 360. ; 360 degrees of longitude correspond to more than one orbit ; because of Earth rotation ylat = inclin * COS( !DTOR * (longt - orblong) / perrat ) subso = SORT(longt) OPLOT, longt(subso), ylat(subso), THICK=2., COLOR=15 ; Find current location mn = MIN( ABS( longt - sclong ), imn ) sclat = ylat(imn) ; Now tailor to exact window size ; Where is center of plot pc = CONVERT_COORD( plong, 0., /DATA, /TO_NORM ) pc = pc(0:1) ; Where should the spacecraft be? sc = CONVERT_COORD( sclong, sclat, /DATA, /TO_NORM ) sc = sc(0:1) ; Increase the offsets because YO_LOGO restricts itself sc(0) = pc(0) + (sc(0) - pc(0)) * ylongscale sc(1) = pc(1) + (sc(1) - pc(1)) * ylatscale YO_LOGO, sc(0), sc(1), yos if data_chk(giffile, /string) then begin zbuff2file, giffile, /gif endif ; BOTTOM OF REPEAT LOOP IF( KEYWORD_SET(replot) EQ 0) THEN RETURN IF( replot EQ 0) THEN RETURN WAIT, 60*replot GOTO, REPEATLOOP END