;$author: $ ;$Date: 2012/05/11 20:54:43 $ ;$Header: /home/cdaweb/dev/control/RCS/plot_enaflux5.pro,v 1.17 2012/05/11 20:54:43 johnson Exp johnson $ ;$Locker: johnson $ ;$Revision: 1.17 $; ; ;------------------------------------------------------------------------- ; 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. ; PNG : Creates PNG file instead of Xwindow. If set to a string, ; this will be the name of the gif, if set to 1, the png ; will be named idl.png. ; 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. ; ;Copyright 1996-2013 United States Government as represented by the ;Administrator of the National Aeronautics and Space Administration. ;All Rights Reserved. ; ;------------------------------------------------------------------ ; 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, PNG=PNG,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 if keyword_set(PNG) then begin s = size(PNG) if (s[n_elements(s)-2] ne 7) then PNG='idl.png' 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 if keyword_set(PNG) then begin xscale=!x.s & yscale=!y.s & bytemap=tvrd() & tvlct,r,g,b,/get s = size(PNG) & ns = n_elements(s) if s[ns-2] eq 7 then gname = PNG else gname = 'idl.png' write_png,gname,bytemap,r,g,b device,/close & set_plot,'X' endif return,0 end PRO testplot_enaflux3,imgfile,GIF=GIF,PNG=PNG,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,PNG=PNG,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,PNG=PNG,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,PNG=PNG,$ 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