;$author: $ ;$Date: 2006-10-11 13:32:51 -0700 (Wed, 11 Oct 2006) $ ;$Header: /home/cdaweb/dev/control/RCS/plot_enaflux5.pro,v 1.13 2006/04/05 14:47:06 klipsch Exp $ ;$Locker: $ ;$Revision: 8 $; ; ;------------------------------------------------------------------------- ; NAME: plot_enaflux ; PURPOSE: ; Plot a given image that is assumed to be square and be centered ; on the Earth. Spacecraft spin axis orientation and orbit data ; are used to scale the Earth and rotate the image so that the ; Earth's north is up. Dipole field lines are also drawn. ; Note that the original input image is scaled (modified) ; by this routine. ; ; CALLING SEQUENCE: ; out = plot_enaflux(etime,image,fov,resolution,sc_pos_gci, ; sc_spinaxis_gci,nadir) ; INPUTS: ; etime = time of the image, in cdf_epoch units ; image = 2d image to be plotted. Must be square image. ; fov = total width of field of view, in degrees. ; resolution = resolution of each pixel in image, in degrees. ; sc_pos_gci = spacecraft position vector [3] in ECI coords. ; sc_spinaxis_gci = spacraft spinaxis uvector [3] in ECI coords. ; nadir = 1=earthpointing, 0=anti-earthpointing ; ; KEYWORD PARAMETERS: ; REFORMER : If image is not perfectly square and this keyword ; is set, the image will be streched to be square. ; ORBIT : In the event that the spacecraft spin axis unit vector ; is not available to the user of this function, and ; given that this plotting code makes the assumption ; that the spin axis is normal to the orbit plane, setting ; this keyword to the [RAAN,INCLINATION] of the orbit ; will cause this routine to compute the sc_spinaxis_gci ; input parameter for the user. ; REVERSEORDER : The 2d image is assumed to be [spinangle (i.e. time) ; by elevation], and that the spinangle and elevation are ; in increasing order. Setting this keyword executes ; the idl reverse,1 which can be used if the spinangle is ; decreasing rather than increasing. Setting this keyword ; to 2 will cause the image to be transposed and reversed. ; GIF : Creates GIF file instead of Xwindow. If set to a string, ; this will be the name of the gif, if set to 1, the gif ; will be named idl.gif. ; WSIZE : Causes the plotcode to apply IDL's CONGRID function ; to change the image size. Example, WSIZE=[200,200] ; NOCIRCLES: If set, equatorial cirlces at 3.3 and 6 Re won't be drawn. ; NODIPOLES: If set, magfield dipoles won't be drawn. ; NOEARTH : If set, the Earth won't be overlaid onto image. ; EDGES : If set to 1, roberts edge enhancement will be applied. ; If set to 2, sobel edge enhancement will be applied. ; NOBORDER : If set, no extra border space will be made around image. ; NOCOLORBAR : If set, no colorbar will be added to window. ; SMOOTH : If set, boxcar average smoothing will be applied. ; SCALEMIN : If set, will be used as the minimum scale instead of 1. ; SCALEMAX : If set, will be used as the maximum scale instead of max(image) ; DEBUG : If set, additional output is printed. ; OUTPUTS: ; out = string array ; AUTHOR: ; Mei-Ching Fok NASA/GSFC Original version December,1999. ; MODIFICATION HISTORY: ; Richard Burley NASA/GSFC/632 plot_enaflux wrapper put around Mei-Chings ; 5/24/2001 original code. Keywords added. ; Tami Kovalick Raytheon ITSS has re-integrated new releases of this software ; into the CDAWeb version of the s/w. ; 4/24/2001 and again 6/19/2001 ; Richard Burley NASA/GSFC/632 Enhanced reverseorder keyword to use ; 6/22/2001 multiple values to apply multiple reversals. ;----------------------------------------------------------------------------- FUNCTION plot_enaflux3,etime,image,spin_angle,polar_angle,np,sc_pos_sm, $ sc_spin_axis_sm,sc_pos_geo,geo_N_sm,nadir,GIF=GIF,WSIZE=WSIZE,$ NOCIRCLES=NOCIRCLES,NODIPOLES=NODIPOLES,NOEARTH=NOEARTH,$ NOCOLORBAR=NOCOLORBAR,EDGES=EDGES,NOBORDER=NOBORDER,$ SMOOTH=SMOOTH,SCALEMAX=SCALEMAX,SCALEMIN=SCALEMIN,$ REVERSEORDER=REVERSEORDER,DEBUG=DEBUG ; Convert time from cdf-epoch to year, day, hour and minute year=0 & month=0 & day=0 & hour=0 & minute=0 & sec=0 cdf_epoch,etime,year,month,day,hour,minute,sec,/break ; Make sure a square image for i=0,np-1 do begin if (spin_angle(i) ne polar_angle(i)) then begin print,' Error: spin angles and polar angles are not the same.' return,-1 endif endfor ; The rotation code in this routine expects the image array order to be ; [spinangle(i.e. time),elevation]. If the image is in the reverse order ; the /REVERSEORDER keyword should be set, and this code will correct. if keyword_set(REVERSEORDER) then begin if REVERSEORDER eq 1 then image = reverse(image,1) if REVERSEORDER eq 2 then image = reverse(transpose(image),2) endif ; Set min/max dimension sizes pix_size=spin_angle(1)-spin_angle(0) x_min=spin_angle(0)-0.5*pix_size x_max=spin_angle(np-1)+0.5*pix_size y_min=x_min & y_max=x_max ; Make sure sc_pos_sm is perpendicular to sc_spin_axis_sm dotproduct=0. for i=0,2 do dotproduct=dotproduct+sc_pos_sm(i)*sc_spin_axis_sm(i) ; Note, RB loosened up the normal requirement from 1.e-4 to 1.e-3 if (abs(dotproduct) gt 1.e-3) then begin print,' Error: spin axis not normal to the spacecraft position vector. Dotp=',dotproduct return,-1 endif ; Find the SM components of x-axis (azimuthal angle) of the satellite frame, in ; which y-axis is along the spin axis and z-axis is along the satellite pos ; vector. For image looking outward from Earth, z-axis is pointed from the ; satellite to the Earth. x_sat=fltarr(3) & y_sat=fltarr(3) & z_sat=fltarr(3) rspin=sqrt(sc_spin_axis_sm(0)*sc_spin_axis_sm(0)+sc_spin_axis_sm(1)* $ sc_spin_axis_sm(1)+sc_spin_axis_sm(2)*sc_spin_axis_sm(2)) rs2=(sc_pos_sm(0)*sc_pos_sm(0)) + (sc_pos_sm(1)*sc_pos_sm(1)) + $ (sc_pos_sm(2)*sc_pos_sm(2)) rs=sqrt(rs2) & rs_tst=sqrt(rs2-1.) z_sat(*)=nadir*sc_pos_sm(*)/rs ; changed to -1 * sc_spin_axis on 2/16 by MC & RB ;y_sat(*)=sc_spin_axis_sm(*)/rspin y_sat(*)= (-1.0 * sc_spin_axis_sm(*))/rspin x_sat(0)=y_sat(1)*z_sat(2)-y_sat(2)*z_sat(1) ; x_sat = y_sat x z_sat x_sat(1)=y_sat(2)*z_sat(0)-y_sat(0)*z_sat(2) x_sat(2)=y_sat(0)*z_sat(1)-y_sat(1)*z_sat(0) ; Rotate satellite frame such that magnetic dipole axis is up (along y-axis) ang_rotate=atan(y_sat(2),x_sat(2))-!pi/2 ; angle between dipole axis & y_sat ;RB commented the following several lines out ;x_sat_new=fltarr(3) & y_sat_new=fltarr(3) & z_sat_new=fltarr(3) ;z_sat_new(*)=z_sat(*) ;x_sat_new(0)=-z_sat(1) ;x_sat_new(1)=z_sat(0) ;x_sat_new(2)=0. ;y_sat_new(0)=z_sat_new(1)*x_sat_new(2)-z_sat_new(2)*x_sat_new(1) ;y_sat_new(1)=z_sat_new(2)*x_sat_new(0)-z_sat_new(0)*x_sat_new(2) ;y_sat_new(2)=z_sat_new(0)*x_sat_new(1)-z_sat_new(1)*x_sat_new(0) ; Rotate and scale the image c_ang_rotate=cos(ang_rotate) & s_ang_rotate=sin(ang_rotate) new_image=fltarr(np,np) ; ; Scale the image - Scaling business added by RB - integrated by TJK on 4/23/01 ; Determine the proper scale for the image flx_min=1.0 & if keyword_set(SCALEMIN) then flx_min=SCALEMIN flx_max=max(image) & if keyword_set(SCALEMAX) then flx_max=SCALEMAX if keyword_set(DEBUG) then print, 'FLX min and max = ',flx_min, flx_max if (flx_max lt flx_min) then begin if keyword_set(DEBUG) then print,'ERROR>plot_enaflux>flx_max < flx_min!' return,-1 endif if (flx_max - flx_min) lt 1.0 then flx_max = flx_max + 1.0 ; make sure enough range log_min=alog10(flx_min) & log_max=alog10(flx_max) if keyword_set(DEBUG) then begin print,'INFO>plot_enaflux>flx_min =',flx_min,' log_min=',log_min print,'INFO>plot_enaflux>flx_max =',flx_max,' log_max=',log_max endif ; THE FOLLOWING 2 LINES WERE COMMENTED OUT ON 2/8 BY RB SO THAT GIFS ; COULD BE AUTOSCALED TO HELP MEI-CHING FIGURE OUT THE ROTATION PROBLEM. ;flx_min=1. & log_min=alog10(flx_min) ; min. flux will be plotted ;flx_max=1.e3 & log_max=alog10(flx_max) ; max. flux will be plotted ; THE FOLLOWING 2 LINES WERE ADDED ON 2/8 BY RB TO AUTOSCALE ;flx_min=1. & log_min=alog10(flx_min) ; min. flux will be plotted ;flx_max=max(image) & log_max=alog10(flx_max) ; max. flux will be plotted ;print,'FLUX_MIN/MAX = ',flx_min,' ',flx_max ;print,'LOG__MIN/MAX = ',log_min,' ',log_max if !d.window ge 0 then loadct, 13, ncolors=240 ;ncolor = !d.table_size - 2 ncolor = 240 - 2 fcolor = float(ncolor) for i=0,np-1 do begin for j=0,np-1 do begin ; --- MCF begin comment out --- ; xn=spin_angle(i)*c_ang_rotate-spin_angle(j)*s_ang_rotate ; xni=(xn-spin_angle(0))/pix_size ; yn=spin_angle(i)*s_ang_rotate+spin_angle(j)*c_ang_rotate ; xnj=(yn-spin_angle(0))/pix_size ; result=interpolate(image,[xni,xni],[xnj,xnj],missing=flx_min) ; if (result(0,0) lt flx_min) then result(0,0)=flx_min ; if (result(0,0) gt flx_max) then result(0,0)=flx_max ; log_flx=alog10(result(0,0)) ; log (flux) ; ; scale the flux to from 1 - fcolor ; new_image(i,j)=1.+(fcolor-1.)*(log_flx-log_min)/(log_max-log_min) value = image[i,j] if value lt flx_min then value = flx_min if value gt flx_max then value = flx_max log_flx = alog10(value) ; scale the flux to from 1 - fcolor image[i,j] = 1.+(fcolor-1.)*(log_flx-log_min)/(log_max-log_min) endfor endfor ; Determine the window size based on keyword or default x_wsize=600 & y_wsize=650 ; original sizes from meiching if keyword_set(WSIZE) then begin ; validate the keyword valid=1 & s=size(WSIZE) & ns=n_elements(s) if (s(0) ne 1) or (s(1) ne 2) then valid = 0 if (s(ns-2) lt 2) or (s(ns-2) gt 3) then valid = 0 if valid eq 1 then begin if (min(wsize) lt 40) or (max(wsize) gt 800) then valid = 0 endif if valid eq 1 then begin x_wsize=wsize(0) & y_wsize=wsize(1) & endif endif if keyword_set(NOBORDER) then begin im_size=fix(x_wsize/np)*np ; make sure im_size is multiple of np x_wsize=im_size & y_wsize=im_size endif ; Plot image to Xwindow or GIF if keyword_set(GIF) then begin s = size(GIF) if (s(n_elements(s)-2) ne 7) then GIF='idl.gif' set_plot,'z' tvlct,red,green,blue,/get ; loadct, 13 ;Rick changed and integrated by TJK 6/19/2001 ; red[0]=255 & green[0]=255 & blue[0]=255 ; make color=0 white ;TJK, change hardcoded array value to !d.n_colors-1 ; red[254]=255 & green[254]=255 & blue[254]=255 & mywhite=254 ; make color=254 white mywhite = !d.table_size-1 red [mywhite] = 255 green[mywhite] = 255 blue [mywhite] = 255 tvlct, red, green, blue device,set_resolution=[x_wsize,y_wsize],set_char=[6,11],z_buffering=0 ;Rick changed and integrated by TJK 6/19/2001 mywhite=0 & myblack=1 ; deviceopen,6,fileOutput=GIF,sizeWindow=[x_wsize,y_wsize] endif else begin ; open x-window display ; loadct, 13 tvlct, red, green, blue, /get ;TJK - not sure about this change... ; red[0]=255 & green[0]=255 & blue[0]=255 ; make color=0 white red[254]=255 & green[254]=255 & blue[254]=255 & mywhite=254 ; make color=254 white tvlct, red, green, blue window,/FREE,xpos=420,ypos=10,xsize=x_wsize,ysize=y_wsize ;TJK not sure about this either mywhite=0 & myblack=1 myblack=0 endelse ; Size and display the image ;This new section was added on 4/23/01 by TJK - this is from RB's latest ;version that he sent to us... x0i=0.1*x_wsize & y0i=0.2*y_wsize ; MCF begin add if keyword_set(NOBORDER) then begin x0i=0. & y0i=0. endif x0=x0i/x_wsize & y0=y0i/y_wsize ; scale image size to 0-1 x1=(x0i+im_size)/x_wsize & y1=(y0i+im_size)/y_wsize ang_rotate_d=ang_rotate*180./!pi ; rotation angle in degree ;print,'INFO>plot_enaflux3>ang_rotate_d=',ang_rotate_d ;print,'ang_rotate_d=',ang_rotate_d ; RBDEBUG map_set,0.,0.,ang_rotate_d,/azimuth,/iso,position=[x0,y0,x1,y1], $ limit=[y_min,x_min,y_max,x_max],/noborder new_image=map_image(image,sx,sy,x_size,y_size,latmin=y_min,$ latmax=y_max,lonmin=x_min,lonmax=x_max) ; MCF end add if not keyword_set(NOBORDER) then im_size=fix(0.8*x_wsize/np)*np c_image=congrid(new_image,im_size,im_size) if keyword_set(EDGES) then begin ; apply edge enhancement if EDGES eq 1 then c_image=roberts(c_image) if EDGES eq 2 then c_image=sobel(c_image) endif ;RB added w/ second version of this s/w ;Rick's version of smoothing that has artifacts... doesn't smooth the whole ;image... if keyword_set(SMOOTH) then begin ; apply boxcar average smoothing if SMOOTH eq 1 then begin ; compute smoothing parameter SMOOTH = ceil(im_size / 7.0) & if (SMOOTH mod 2) eq 0 then SMOOTH=SMOOTH+1 endif print, 'Images being SMOOTHED by a factor of, ',smooth c_image = smooth(c_image,smooth,/edge_truncate) ; smooth factor changes based on the ; image size endif ;TJK my version follows, which doesn't have artifacts (no unsmoothed blocks ;around the edges. ;if keyword_set(SMOOTH) then begin ; apply boxcar average smoothing ; if SMOOTH eq 1 then begin ; compute smoothing parameter ; SMOOTH = ceil(im_size / 10.0) & if (SMOOTH mod 2) eq 0 then SMOOTH=SMOOTH+1 ; endif ; c_image = smooth(c_image,10) ;changed to a hard number ; print, 'Images being SMOOTHED by a factor of 10' ;endif ;MCF commented out in second version ;x0i=0.1*x_wsize & y0i=0.2*y_wsize ;if keyword_set(NOBORDER) then begin ; x0i=0. & y0i=0. ;endif ;x0=x0i/x_wsize & y0=y0i/y_wsize ; scale image size to 0-1 ;x1=(x0i+im_size)/x_wsize & y1=(y0i+im_size)/y_wsize ;npt=480 ; no. of points in draw fieldline, circle.. ;End of MCF comment out npt=1000 ; no. of points in draw fieldline, circle.. ;MCF add e_rad = asin(1.0/rs) * 180.0 / !pi ; angle subtaned by the Earth e_x = fltarr(npt+1) & e_y = fltarr(npt+1) x = fltarr(2) & y = fltarr(2) ;MCF changed tv,c_image,x0i,y0i ; display the image tv,c_image,sx/200,sy/200 ; display the image ; MCF add ; RB, added division ; by 200 for Z-buffer ; Add color bar and label if not keyword_set(NOCOLORBAR) and not keyword_set(NOBORDER) then begin tempvar = bytarr(ncolor,2) for i=0,ncolor-1 do begin tempvar(i,0) = i+1 & tempvar(i,1) = i+1 endfor color_bar=congrid(tempvar,ncolor,30) tv,color_bar,x0i,0.32*y0i xyouts,x0,0.2*y0,string(log_min,'(f3.1)'),$ alignment=0.5,size=1.5,/normal,color=myblack xlab = x0+float(ncolor)/x_wsize xyouts,xlab,0.2*y0,string(log_max,'(f3.1)'),$ alignment=0.5,size=1.5,/normal,color=myblack xyouts,0.5*(xlab+x0),0.1*y0,'log (particle/cm2/sr/s)',$ alignment=0.5,size=1.5,/normal,color=myblack endif ;x1=(x0i+im_size)/x_wsize ;y1=(y0i+im_size)/y_wsize ;npt=480 ; no. of points in draw fieldline, circle.. ;e_rad = asin(1.0/rs) * 180.0 / !pi ; angle subtaned by the Earth ;e_x = fltarr(npt+1) & e_y = fltarr(npt+1) ;x = fltarr(2) & y = fltarr(2) ; Calculate and draw circles at 3 and 6.6 Re at the equator and ; draw connections between them r=fltarr(3) & rp=r ; MCF add three_hr=!pi/4. & ang0=atan(sc_pos_sm(1),sc_pos_sm(0))+!pi del = 2.0 * !pi / npt & ncpt=16 xc=fltarr(ncpt) & yc=fltarr(ncpt) & lt_lab=strarr(ncpt) for i=1,2 do begin & localtime0=0 ; init if (i eq 1) then rc=3. & if (i eq 2) then rc=6.6 icn=i-1 & ic=-1 for ii = 0,npt do begin & ibehind=nadir ; init ang = ang0+ii*del ; xe = rc * cos(ang) ;MCF begin comment out ; ye = rc * sin(ang) & r2=xe*xe+ye*ye ; re2=(sc_pos_sm(0)-xe)^2+(sc_pos_sm(1)-ye)^2+sc_pos_sm(2)^2 ; re=sqrt(re2) ; dist. ./. satellite. & equator ; cosa=nadir*(rs2+re2-r2)/(2.*rs*re) ; angl=acos(cosa)*180./!pi ;angle subtended in deg ; xd=xe*x_sat_new(0)+ye*x_sat_new(1) ; yd=xe*y_sat_new(0)+ye*y_sat_new(1) ; rd=sqrt(xd*xd+yd*yd) ;MCF end comment out r(0) = rc * cos(ang) ; MCF begin add r(1) = rc * sin(ang) r(2) = 0. rp(*)=r(*)-sc_pos_sm(*) re=sqrt(total(rp*rp)) ; dist. ./. satellite. & equator xp=total(rp*x_sat) yp=total(rp*y_sat) zp=total(rp*z_sat) xd=atan(xp,-zp)*180./!pi yd=asin(yp/re)*180./!pi angl=sqrt(xd*xd+yd*yd) ;angle subtended in deg ; MCF edd add if (angl gt e_rad or re lt rs_tst or nadir eq -1) then begin ibehind=-1 ; this point can be seen ; ic=ic+1 & e_x(ic)=xd*angl/rd & e_y(ic)=yd*angl/rd ;MCF comment out ic=ic+1 & e_x(ic)=xd & e_y(ic)=yd ; MCF add endif ang_o_3=fix(ang/three_hr) & test=three_hr*ang_o_3 if ((ang-test) lt del and icn le (ncpt-1)) then begin ; 3hr Label localtime=ang_o_3*3-12 ; -12 cause 180 deg -> 00 LT if (localtime lt 0) then localtime=localtime+24 if (localtime gt 24) then localtime=localtime-24 if (icn eq (i-1) or localtime ne localtime0) then begin ;avoid reps lt_lab(icn)=string(localtime,'(i2.2)') if (ibehind eq -1) then begin xc(icn)=e_x(ic) & yc(icn)=e_y(ic) endif else begin ; if behind earth, put point on earth rim ; xc(icn)=xd*e_rad/rd & yc(icn)=yd*e_rad/rd replaced w/ following ; line on 12/28 based on email from mei-ching. xc(icn)=xd*e_rad/angl & yc(icn)=yd*e_rad/angl endelse icn=icn+2 endif localtime0=localtime endif endfor ;TJK, change color keyword value below to !d.n_colors-1 instead of zero ;change all "color=0" to color=!d.n_colors-1 ;change all settings of xticks=8,yticks=8 to xticks=1,yticks=1 to suppress ;tick marks. Change the x/ystyle to +4 so that the axis won't be drawn. ;MCF commented this section out and replaced it w/ one line... ; if (i eq 1) then begin ; if not keyword_set(NOCIRCLES) then begin ; plot,e_x(0:ic),e_y(0:ic),position=[x0,y0,x1,y1],$ ; xrange=[x_min,x_max],yrange=[y_min,y_max],xticks=1,yticks=1,$ ; xstyle=1+4,ystyle=1+4,color=!d.n_colors-1,/noerase ; endif ; endif ; if (i eq 2) then if not keyword_set(NOCIRCLES) then $ ; oplot,e_x(0:ic),e_y(0:ic),color=!d.n_colors-1 ;TJK - mywhite is set to 0 above, MCF had color set to mywhite - I would prefer ; the use of !d.n_colors-1 instead ;if not keyword_set(NOCIRCLES) then oplot,e_x(0:ic),e_y(0:ic),color=mywhite ; MCF add if not keyword_set(NOCIRCLES) then $ oplot, e_x(0:ic), e_y(0:ic), color=!d.table_size-1 endfor for i=0,14,2 do begin if not keyword_set(NOCIRCLES) then $ oplot, xc(i:i+1), yc(i:i+1), color=!d.table_size-1 labx=1.1*xc(i+1) & laby=1.1*yc(i+1) if (abs(labx) le x_max and abs(laby) le x_max) then begin if not keyword_set(NOCIRCLES) then $ xyouts, labx, laby, lt_lab(i+1), color=!d.table_size-1 endif endfor ; Calculate and draw dipole fieldlines at L=3, 6.6, MLT=0, 6, 12, 18 for i=1,2 do begin if (i eq 1) then rc=3. & if (i eq 2) then rc=6.6 sinang0=sqrt(1./rc) & ang0=asin(sinang0) ; draw from earth surface del=(!pi-2.*ang0)/npt for j=0,270,90 do begin phi=j*!pi/180. ; azimuthal angle in radian cphi=cos(phi) & sphi=sin(phi) for n=0,1 do begin ; do 2 hemispheres separately i1=n*npt/2 & i2=i1+npt/2 & ic=-1 for ii=i1,i2 do begin ang=ang0+ii*del & sang=sin(ang) ;MCF commented out ; r=rc*sang*sang & r2=r*r ; xe=r*sang*cphi & ye=r*sang*sphi & ze=r*cos(ang) ; re2=(sc_pos_sm(0)-xe)^2+(sc_pos_sm(1)-ye)^2+(sc_pos_sm(2)-ze)^2 ; re=sqrt(re2) ;dist. ./. sat. & fieldline ; cosa=nadir*(rs2+re2-r2)/(2.*rs*re) ; angl=acos(cosa)*180./!pi ;angle subtended in deg ; xd=xe*x_sat_new(0)+ye*x_sat_new(1)+ze*x_sat_new(2) ; yd=xe*y_sat_new(0)+ye*y_sat_new(1)+ze*y_sat_new(2) ; rd=sqrt(xd*xd+yd*yd) r1=rc*sang*sang ; dipole fieldline equation ; MCF begin add r2=r1*r1 r(0)=r1*sang*cphi r(1)=r1*sang*sphi r(2)=r1*cos(ang) rp(*)=r(*)-sc_pos_sm(*) re=sqrt(total(rp*rp)) ; dist. ./. satellite.&fieldlines xp=total(rp*x_sat) yp=total(rp*y_sat) zp=total(rp*z_sat) xd=atan(xp,-zp)*180./!pi yd=asin(yp/re)*180./!pi angl=sqrt(xd*xd+yd*yd) ;angle subtended in deg ; MCF end add if (angl gt e_rad or re lt rs_tst or nadir eq -1) then begin ic=ic+1 ; fieldline point can be seen ; e_x(ic)=xd*angl/rd & e_y(ic)=yd*angl/rd ;MCF comment out e_x(ic)=xd & e_y(ic)=yd ; MCF add endif endfor ;TJK, change color keyword value below to !d.n_colors-1 instead of zero ;change all "color=0" to color=!d.n_colors-1 ;change all settings of xticks=8,yticks=8 to xticks=1,yticks=1 to suppress ;tick marks. Change the x/ystyle to +4 so that the axis won't be drawn. if not keyword_set(NODIPOLES) then begin if (ic gt 0) then begin if keyword_set(NOCIRCLES) then begin ; need to call plot plot,e_x(0:ic),e_y(0:ic),position=[x0,y0,x1,y1],$ xrange=[x_min,x_max],yrange=[y_min,y_max],xticks=1+4,$ yticks=1+4,xstyle=1,ystyle=1,$ color=!d.table_size-1,/noerase endif else oplot, e_x(0:ic),e_y(0:ic), color=!d.table_size-1 endif endif endfor endfor endfor ; Add direction and label of the dipole axis when nadir = 1 if (nadir eq 1) then begin x[0] = 0. & y[0] = 0. & x[1] = x[0] & y[1] = 0.8*x_max ;oplot, x,y, color=0 ;xyouts,x[1],y[1], 'MagDipole', color=0 endif ;TJK, change all settings of xticks=8,yticks=8 to xticks=1,yticks=1 to ;suppress tick marks. ;change color keyword value below to !d.n_colors-1 instead of zero ;change all "color=0" to color=!d.n_colors-1 ;Change the x/ystyle to +4 so that the axis won't be drawn. ; Label the image if (keyword_set(NOCIRCLES)) and (keyword_set(NODIPOLES)) then begin ; Need to call plot with /nodata to set the plot scale since it hasn't been done yet. plot,[x_min,x_max],[y_min,y_max],/nodata,xstyle=1+4,ystyle=1+4,$ xticks=1,yticks=1,color=!d.table_size-1,/noerase,$ position=[x0,y0,x1,y1] endif ;xy outs were commented out on 2/21 by RB ;dang=x_max/2 & base=-1.055*x_max ;for i=-2,2 do begin ; angle=i*dang ; xyouts,angle,base,string(abs(angle),'(i2)'),alignment=0.5,color=1 ; xyouts,base,angle,string(abs(angle),'(i2)'),alignment=0.5,color=1 ;endfor ;xyouts,0.,1.07*base,'degree',alignment=0.5,color=1 ;xyouts,1.07*base,-0.07*x_max,'degree',orientation=90,color=1 ;TJK, change color keyword value below to !d.n_colors-1 instead of zero ;change all "color=0" to color=!d.n_colors-1 ;change all settings of xticks=8,yticks=8 to xticks=1,yticks=1 to suppress ;tick marks. Change the x/ystyle to +4 so that the axis won't be drawn. ; Add Earth outline and continents when nadir = 1 del = 2.0 * !pi / npt if (nadir eq 1) then begin for ii = 0,npt do begin ang = float(ii) * del e_x(ii) = e_rad * cos(ang) & e_y(ii) = e_rad * sin(ang) endfor if not keyword_set(NOEARTH) then begin if keyword_set(NOCIRCLES) then begin ; must call plot instead of oplot plot,e_x,e_y,position=[x0,y0,x1,y1],$ xrange=[x_min,x_max],yrange=[y_min,y_max],xticks=1,yticks=1,$ xstyle=1+4,ystyle=1+4,color=!d.table_size-1,/noerase endif else oplot,e_x,e_y,color=!d.table_size-1 endif ;MCF commented out ; xm=(x1+x0)/2. & e_radx=e_rad*(x1-x0)/(x_max-x_min) ; ym=(y1+y0)/2. & e_rady=e_rad*(y1-y0)/(y_max-y_min) xm=(x1+x0)/2. & e_radx=e_rad*(x1-x0)*2/(180.+x_max-x_min) ; MCF begin add ym=(y1+y0)/2. & e_rady=e_rad*(y1-y0)*2/(180.+y_max-y_min) ; MCF end add geo_lat=asin(sc_pos_geo(2)/rs)*180./!pi ; geographic lat (deg) geo_lon=atan(sc_pos_geo(1),sc_pos_geo(0))*180./!pi ; geographic lon (deg) ;MCF commented out ; xN=geo_N_sm(0)*x_sat_new(0)+geo_N_sm(1)*x_sat_new(1)+geo_N_sm(2)*x_sat_new(2) ; yN=geo_N_sm(0)*y_sat_new(0)+geo_N_sm(1)*y_sat_new(1)+geo_N_sm(2)*y_sat_new(2) ; gamma=90. - atan(yN,xN)*180./!pi ; rotation in deg, clockwise from north xN=geo_N_sm(0)*x_sat(0)+geo_N_sm(1)*x_sat(1)+geo_N_sm(2)*x_sat(2) ;MCF begin add yN=geo_N_sm(0)*y_sat(0)+geo_N_sm(1)*y_sat(1)+geo_N_sm(2)*y_sat(2) gamma=ang_rotate_d+90.-atan(yN,xN)*180./!pi ; rotation(deg), clockwise from N ;MCF end add if not keyword_set(NOEARTH) then begin map_set,geo_lat,geo_lon,$ position=[xm-e_radx,ym-e_rady,xm+e_radx,ym+e_rady],$ /satellite,sat_p=[rs,0.,gamma],/continents,$ con_color=!d.table_size-1,/noborder,/noerase endif endif ; Close the GIF file if writing to GIF if keyword_set(GIF) then begin xscale=!x.s & yscale=!y.s & bytemap=tvrd() & tvlct,r,g,b,/get s = size(GIF) & ns = n_elements(s) if s(ns-2) eq 7 then gname = GIF else gname = 'idl.gif' write_gif,gname,bytemap,r,g,b device,/close & set_plot,'X' endif return,0 end PRO testplot_enaflux3,imgfile,GIF=GIF,WSIZE=WSIZE,NOCIRCLES=NOCIRCLES,$ NODIPOLES=NODIPOLES,NOEARTH=NOEARTH,$ NOCOLORBAR=NOCOLORBAR,EDGES=EDGES,NOBORDER=NOBORDER m=0 & d= 0 & monday,2000,100,m,d cdf_epoch,etime,2000,m,d,/compute & etime=etime+48600000 np=32 & image=fltarr(np,np) spin_angle=intarr(np) & polar_angle=intarr(np) for i=0,31 do begin ; compute centers of 4x4 pixels spin_angle(i)=-62+i*4 & polar_angle(i)=spin_angle(i) endfor nadir=1 ; 1=earthward view, 0=anti-earthward view openr,1,imgfile sc_pos_sm=fltarr(3) & s='' & readf,1,s & reads,s,sc_pos_sm sc_spin_axis_sm=fltarr(3) & s='' & readf,1,s & reads,s,sc_spin_axis_sm sc_pos_geo=fltarr(3) & s='' & readf,1,s & reads,s,sc_pos_geo geo_N_sm=fltarr(3) & s='' & readf,1,s & reads,s,geo_N_sm np=32 & image=fltarr(np,np) & s=strarr(205) & readf,1,s & close,1 s2='' & for i=0,204 do s2=s2+s(i) & reads,s2,image s=plot_enaflux3(etime,image,spin_angle,polar_angle,np,sc_pos_sm,$ sc_spin_axis_sm,sc_pos_geo,geo_N_sm,nadir,$ GIF=GIF,WSIZE=WSIZE,NOCIRCLES=NOCIRCLES,NODIPOLES=NODIPOLES,$ NOEARTH=NOEARTH,NOCOLORBAR=NOCOLORBAR,EDGES=EDGES,$ NOBORDER=NOBORDER) end FUNCTION plot_enaflux,etime,image,fov,resolution,sc_pos_gci,sc_spinaxis_gci,$ nadir,REFORMER=REFORMER,ORBIT=ORBIT,GIF=GIF,WSIZE=WSIZE,$ NOCIRCLES=NOCIRCLES,NODIPOLES=NODIPOLES,NOEARTH=NOEARTH,$ NOCOLORBAR=NOCOLORBAR,EDGES=EDGES,NOBORDER=NOBORDER,$ SMOOTH=SMOOTH,SCALEMIN=SCALEMIN,SCALEMAX=SCALEMAX,$ REVERSEORDER=REVERSEORDER,DEBUG=DEBUG ; Convert spacecraft position vector (GCI) and spinaxis (GCI) to SM coords. year=0 & month=0 & day=0 & hour=0 & minute=0 & sec=0 ; init params for recalc recalc,year,day,hour,min,sec,epoch=etime ; setup conversion values ; Create scalar variables required when calling geopack routines xgci=sc_pos_gci(0) & ygci=sc_pos_gci(1) & zgci=sc_pos_gci(2) xgse=0.0 & ygse=0.0 & zgse=0.0 & xgsm=0.0 & ygsm=0.0 & zgsm=0.0 xsm=0.0 & ysm=0.0 & zsm=0.0 & xgeo=0.0 & ygeo = 0.0 & zgeo=0.0 ; ; Perform conversions geigse,xgci,ygci,zgci,xgse,ygse,zgse,1,etime gsmgse,xgsm,ygsm,zgsm,xgse,ygse,zgse,-1 smgsm,xsm,ysm,zsm,xgsm,ygsm,zgsm,-1 geigeo,xgci,ygci,zgci,xgeo,ygeo,zgeo,1,epoch=etime sc_pos_sm = [xsm,ysm,zsm] / 6378.14 ; convert to Re sc_pos_geo = [xgeo,ygeo,zgeo] / 6378.14 ; convert to Re ; In the event that the spacecraft spin axis unit vector is was not ; available to the user of this function, and given that the plotting code ; makes the assumption that the spin axis is normal to the orbit plane, ; compute what spin_axis should be given the orbit described with right ; ascension of the ascending node and inclination. if keyword_set(ORBIT) then begin raan = orbit(0) * !dtor & incl = orbit(1) * !dtor target_ra = raan - (90.0 * !dtor) & target_dec = (90.0 * !dtor) - incl x = cos(target_dec) * cos(target_ra) y = cos(target_dec) * sin(target_ra) z = sin(target_dec) m = (x^2 + y^2 + Z^2)^0.5 sc_spinaxis_gci = [x/m,y/m,z/m] * !radeg endif ; Convert spacecraft spin axis unit vector from gci to sm. xgci=sc_spinaxis_gci(0) & ygci=sc_spinaxis_gci(1) & zgci=sc_spinaxis_gci(2) xgse=0.0 & ygse=0.0 & zgse=0.0 & xgsm=0.0 & ygsm=0.0 & zgsm=0.0 xsm=0.0 & ysm=0.0 & zsm=0.0 ; Convert spacecraft position vector from gci to sm geigse,xgci,ygci,zgci,xgse,ygse,zgse,1,etime gsmgse,xgsm,ygsm,zgsm,xgse,ygse,zgse,-1 smgsm,xsm,ysm,zsm,xgsm,ygsm,zgsm,-1 ; Convert units from kilometers to earth_radii re = 6378.14 & sc_spinaxis_sm = [xsm/re,ysm/re,zsm/re] ; if producing debug output then output a dot product analysis if keyword_set(DEBUG) then begin dot=fltarr(3) & for i=0,2 do dot(i) = sc_spinaxis_sm(i) * sc_pos_sm(i) print,'INFO>plot_enaflux>DOT(pos,spin)=',dot,' = ',total(dot) endif ; Compute geographic northpole in solar magnetospheric coordinates geogsm,0.0,0.0,1.0,xgsm,ygsm,zgsm,1 ; convert northpole to gsm smgsm,xsm,ysm,zsm,xgsm,ygsm,zgsm,-1 ; convert to sm geo_N_sm = [xsm,ysm,zsm] / re ; Determine the size of the image and verify that it is square s = size(image) & ns = n_elements(s) if s(0) ne 2 then begin print,'ERROR>image parameter is not two dimensional' & return,-1 endif if s(1) ne s(2) then begin if not keyword_set(REFORMER) then begin print,'ERROR>image parameter is not square and /REFORMER keyword not set' return,-1 endif else begin np = max([s(1),s(2)]) & image = reform(image,np,np) endelse endif else np = s(1) ; Given the field of view and angular resolution (both degrees) of the ; instrument which took the image, and assuming that the image is centered ; on the earth, derive the viewing angles in X (spin) and Y (polar). ;spin_angle = intarr(np) & polar_angle = intarr(np) ;start_spin_angle = -1 * fix(0.5 * fov) ;for i=0,np-1 do begin ; spin_angle(i) = start_spin_angle + (i*resolution) + (0.5 * resolution) ; polar_angle(i) = spin_angle(i) ; valid estimate because of square image ;endfor spin_angle = fltarr(np) & polar_angle = fltarr(np) start_spin_angle = -1. * fix(0.5 * fov) for i=0,np-1 do begin spin_angle(i) = start_spin_angle + (i * resolution) + (0.5 * resolution) polar_angle(i) = spin_angle(i) ; valid estimate because of square image endfor ; ; Generate debug output for MeiChing if keyword_set(DEBUG) then begin print,'INFO>plot_enaflux3>sc_pos_sm=',sc_pos_sm print,'INFO>plot_enaflux3>sc_spinaxis_sm=',sc_spinaxis_sm print,'INFO>plot_enaflux3>sc_pos_geo=',sc_pos_geo print,'INFO>plot_enaflux3>geo_N_sm=',geo_N_sm endif ; ; Generate the plot s = plot_enaflux3(etime,image,spin_angle,polar_angle,np,sc_pos_sm,$ sc_spinaxis_sm,sc_pos_geo,geo_N_sm,nadir,GIF=GIF,$ WSIZE=WSIZE,NOCIRCLES=NOCIRCLES,NODIPOLES=NODIPOLES,$ NOEARTH=NOEARTH,NOCOLORBAR=NOCOLORBAR,EDGES=EDGES,$ NOBORDER=NOBORDER,SMOOTH=SMOOTH,DEBUG=DEBUG,$ SCALEMAX=SCALEMAX,SCALEMIN=SCALEMIN,REVERSEORDER=REVERSEORDER) return,s end