;+ ; NAME: TMAGINT ; ; PURPOSE: ; To display images that have their backgrounds removed (by ; using a minimum image) and the intensity of the images ; within 1 deg of the magnetic zenith. ; CATEGORY: ; ; CALLING SEQUENCE: ; TMAGINT,YEAR,MONTH,DAY,CN,FN,HH,MM,SS,SPAN $ ; [,WAIT=WAIT] [,EXP=EXP] [,CAM=CAM] [,FILT=FILT] $ ; [,VRHIJN=VRHIJN] [,COL_TABLE=COL_TABLE] [,double=double] ; INPUTS: ; YEAR,MONTH,DAY: The year month and day that the images were ; produced. ; CN,FN: The camera and filter number. ; HH,MM,SS: The start time, in UT, of the interval from which ; images will be displayed. ; SPAN: The length, in hours, from the start time that images ; will be displayed. ; KEYWORD PARAMETERS: ; /WAIT: Allows user to manually display the images. By default, ; the procedure will cycle through the images as fast as ; possible. ; EXP: The exposure time of the images. By default, 60 seconds ; will be chosen. ; CAM: The camera number of the minimum image that was used ; for subtraction. If not set, then by default CN will ; be used. ; FILT: The filter number of the minimum image that was used ; for subtraction. If not set, then by default FN will ; be used. ; /VRHIJN: Use van Rhijn layers as the minimum images for background ; subtraction. If not set, then the background subtraction ; images will be the normal minimum images. ; COL_TABLE: Color table to be used in displaying and printing ; images. Default is 13. ; /DOUBLE: Displays another set of images. User will be prompted ; to enter information about the second image set. ; OUTPUTS: ; Displays images to WINDOW 0. ; COMMON BLOCKS: ; None. ; SIDE EFFECTS: ; User can, interactively, do one or more of the following: ; - create postscript files of any image; ; - view and plot a meridian slice across an image or ; pair of images at a specified azimuth; ; - extract the mean intensity over a 1x1 deg field at ; any point in an image or pair of images. ; RESTRICTIONS: ; None known. ; PROCEDURE: ; ... ; MODIFICATION HISTORY: ; Written August 1994, by T A Oliynyk. ; Modified extensively October 1994, by D P Steele. ;- PRO tmagint,year,month,day,cn,fn,hh,mm,ss,span,wait=wait,col_table=col_table,$ exp=exp,cam=cam,filt=filt,vrhijn=vrhijn,double=double IF N_PARAMS() EQ 0 THEN BEGIN doc_library,'tmagint' RETURN ENDIF IF N_PARAMS() LT 9 THEN list=1 !ORDER=0 IF NOT KEYWORD_SET(exp) THEN exp=60 sy=STRING(year MOD 100,FORMAT='(I2.2)') sm=STRING(month,FORMAT='(I2.2)') sd=STRING(day,FORMAT='(I2.2)') date=sy+'/'+sm+'/'+sd scn=STRING(cn,FORMAT='(I1)') sfn=STRING(fn,FORMAT='(I1)') path='/jasper/cnsr3_data1/save/' filename=path+sy+sm+sd+scn+sfn+'.bgs' ; find the appropriate .bgs (primary) images tselect_bgsim,keep,utsec,year,month,day,cn,fn,hh,mm,ss,span,$ list=list,exp=exp,cam=cam,filt=filt,vrhijn=vrhijn IF keep(0) EQ -1 THEN RETURN OPENR,unit,filename,/GET_LUN files=ASSOC(unit,BYTARR(66048)) ; find the pixels that are located within 1 deg of the magnetic zenith deg=1 tmagzen,year,month,day,cn,fn,deg,mask,cp,rp,/small,/height index=WHERE(mask EQ 1) ; index to the pixels within 1 deg of the magnetic zenith ; locate the bounds of the region (1 deg of the magnetic zenith) rows=index/256 cols=index MOD 256 minr=MIN(rows,MAX=maxr) minc=MIN(cols,MAX=maxc) ; convert row and column numbers to azimuth and zenith angles column=INDGEN(256)#REPLICATE(1,1,256) row=TRANSPOSE(column) cr2az,year,month,day,cn,fn,column,row,az,za,/small,/deg,site=3 ; get user to decide on a second set of images to display IF KEYWORD_SET(double) THEN BEGIN col_table2=0 decis='' PRINT,"Select camera and filter number of companion images." READ,cn2,fn2 PRINT," " PRINT,"Select the color table that you wish to view the images in." READ,col_table2 PRINT,' ' PRINT,"Do you want a list of available image sets? ([y]/n)" READ,decis decis=STRUPCASE(decis) PRINT,' ' IF decis EQ 'N' THEN BEGIN list=0 PRINT,"Enter the start time (hour,minute,second) of image set." READ,hh2,mm2,ss2 PRINT,"Enter the time (in hours) that the image set spans." READ,span2 PRINT,' ' PRINT,"Do you want the minimum images used for background subtraction to be from" PRINT,"the same camera/filter combination as the unsubtracted images? ([y]/n)" READ,decis decis=STRUPCASE(decis) IF decis EQ 'N' THEN BEGIN PRINT, "Enter the camera and filter number of the minimum image." READ,cam2,filt2 ENDIF ELSE BEGIN cam2=cn2 filt2=fn2 ENDELSE ENDIF ELSE list=1 scn2=STRING(cn2,FORMAT='(I1)') sfn2=STRING(fn2,FORMAT='(I1)') filename2=path+sy+sm+sd+scn2+sfn2+'.bgs' tselect_bgsim,keep2,utsec2,year,month,day,cn2,fn2,hh2,mm2,ss2,span2,$ list=list,exp=exp,cam=cam2,filt=filt2,vrhijn=vrhijn IF keep(0) EQ -1 THEN BEGIN FREE_LUN,unit RETURN ENDIF OPENR,unit2,filename2,/GET_LUN files2=ASSOC(unit2,BYTARR(66048)) ; find the pixels that are containted within 1 deg of the magnetic zenith ; on the secondary images deg2=1. tmagzen,year,month,day,cn2,fn2,deg2,mask2,cp2,rp2,/small,/height index2=WHERE(mask2 EQ deg2) ; index to the pixels within 1 deg of the magnetic zenith ; locate the bounds of the region ( 1 deg of the magetic zenith) on the secondary images rows2=index2/256 cols2=index2 MOD 256 minr2=MIN(rows2,MAX=maxr2) minc2=MIN(cols2,MAX=maxc2) ; convert the row and column numbers to azimuth and elevation angles ; (for secondary images) cr2az,year,month,day,cn2,fn2,column,row,az2,za2,/small,/deg,site=3 ; convert the zenith angles to latitude angles (in radians) gamma=(90-za2)*!DTOR alpha=SIN(gamma) ; compute the SIN of the latitude angles beta=COS(gamma) ; compute the COS of the latitude angles ; if two images are to be displayed, then make window the appropriate size WINDOW,XSIZE=1070,YSIZE=700 xsh=710 & ysh=22 ysh2=44+317 ENDIF ELSE BEGIN ; if only one image is to be displayed, then make the window the appropriate size WINDOW,XSIZE=1070 xsh=710 ysh=185 decis='' ENDELSE ; load color table 13 as a default IF NOT KEYWORD_SET(col_table) THEN col_table=13 IF KEYWORD_SET(double) THEN BEGIN LOADCT,col_table,/SILENT tvlct,r1,g1,b1,/get LOADCT,col_table2,/SILENT tvlct,r2,g2,b2,/get nsteps=N_ELEMENTS(r1) r=REBIN([r1,r2],nsteps,/sample) g=REBIN([g1,g2],nsteps,/sample) b=REBIN([b1,b2],nsteps,/sample) tvlct,r,g,b ENDIF ELSE LOADCT,col_table,/silent ; if two images are going to be displayed, only take images that were produced ; at the same time IF KEYWORD_SET(double) THEN BEGIN sametime=where_array(utsec2(keep2),utsec(keep)) sametime2=where_array(utsec(keep),utsec2(keep2)) IF (sametime(0) EQ -1) OR (sametime2(0) EQ -1) THEN BEGIN PRINT,"None of the images from the two sets were produced at the same time." FREE_LUN,unit FREE_LUN,unit2 RETURN ENDIF sametime=sametime(uniq(sametime,SORT(sametime))) sametime2=sametime2(uniq(sametime2,SORT(sametime))) keep=keep(sametime) ; index into primary images that occur at the same time ; as the secondary images keep2=keep2(sametime2) ; index into secondary images that occur at the same time ; as the primary images ENDIF n=N_ELEMENTS(keep) ; number of images to display ps_plot=[-1] ; initialize arrays to hold information about the primary image del=FLTARR(n) ; array to hold the deviation (in degrees) from the magnetic zenith zen=del ; ... zenith angle of chosen region (either 1 deg about magnetic zenith ; or a 5 by 5 pixel area about a position chosen by the user) azm=del ; ... azimuth angle of region ... var=del ; the standard deviation in the intensity in the region u=del ; the average intensity in the region ccol=INTARR(n); array to hold the column pixels of the center of the region crow=ccol ; array to hold the row pixels of the center region rdim=BYTARR(n); array containing 0's and 1's. 0 if the region is centered on the ; magnetic zenith, and 1 if it is centered somewhere else level=del ; array to hold the threshold level ; zenith and azimuth angle, in degrees, of the magnetic zenith zenith=za(cp,rp) azimuth=az(cp,rp) ; convert the magnetic zenith and azimuth angles to east longitude and latitude ; angles in radians ( for the primary image) lam2=(90-zenith)*!DTOR phi2=azimuth*!DTOR ; rotate the arrays containing zenith and azimuth angles (for the primary image), ; because the image will be rotated later before displaying (to put north facing upwards). az=ROTATE(az,2) za=ROTATE(za,2) IF KEYWORD_SET(double) THEN BEGIN ; initialize arrays to hold information about the secondary image del2=FLTARR(n) ; array to hold the deviation (in degrees) from the magnetic zenith zen2=del2 ; ... zenith angle of chosen region ( either 1 deg about magnetic zenith ; or a 5 by 5 pixel area about a position chosen by the user azm2=del2 ; ... azimuth angle of region ... var2=del2 ; the standard deviation in the intensity in the region u2=del2 ; the average intensity in the region ccol2=INTARR(n) ; array to hold the column pixels of the center of the region crow2=ccol2 ; array to hold the row pixels of the center region rdim2=BYTARR(n) ; array containing 0's and 1's. 0 if the region is centered on the ; magnetic zenith, and 1 if it is centered somewhere else ; zenith and azimuth angle, in degrees, of the magnetic zenith ; find the azimuth and zenith angle of the magnetic zenith for the secondary image zenith2=za2(cp2,rp2) azimuth2=az2(cp2,rp2) ; convert the magnetic zenith and azimuth angles to east longitude and latitude ; angles in radians. lam22=(90-zenith2)*!DTOR phi22=azimuth2*!DTOR az2=ROTATE(az2,2) za2=ROTATE(za2,2) gamma=(90-za2)*!DTOR alpha=SIN(gamma) ; compute the SIN of the latitude angles beta=COS(gamma) ; compute the COS of the latitude angles ENDIF hold=200 csize=1.2 ; set the character size on the display IF KEYWORD_SET(wait) THEN PRINT, 'Type ? for help' i=-1 chgthresh=0 ; flag to determine if a threshold level is set norefresh=0B ; flag suppressing display refresh if set k='' !ORDER=0 REPEAT BEGIN i=FIX((i+1) MOD n) j=keep(i) IF KEYWORD_SET(double) THEN j2=keep2(i) IF NOT norefresh THEN BEGIN rdim(i)=0B del(i)=0. zen(i)=zenith azm(i)=azimuth min_c=minc max_c=maxc min_r=minr max_r=maxr ; find the time the primary images were produced hour=utsec(j)/3600L minute=(utsec(j)-hour*3600L)/60L second=utsec(j)-hour*3600L-minute*60L time=STRING(hour,minute,second,FORMAT='(I2.2,":",I2.2,":",I2.2)') im=BYTE(files(j),512,256,256) ; if the primary images were normalized to 8 bits, then reverse this hdb=BYTE(files(j),0,512) hd=gethd(hdb) power=hd.proc.base2 scale=hd.proc.scale_factor IF scale GT 0 THEN BEGIN IF power EQ 1 THEN im=LONG(im)*(2L^scale) IF power EQ 0 THEN im=LONG(im)*scale ENDIF IF chgthresh EQ 1 THEN BEGIN im=im(-1)) 'B': i=FIX((i-2)>(-1)) 'P': i=FIX((i-2)>(-1)) '+': i=FIX(i) ' ': i=FIX(i) 'I': BEGIN PRINT,STRTRIM(i,2) i=FIX(i-1) END 'N': BEGIN PRINT,STRTRIM(n,2) i=FIX(i-1) END 'F': BEGIN wait=0B i=FIX(i-1) END 'S': BEGIN wait=1B i=FIX(i-1) END 'T': BEGIN thresh=0 PRINT,FORMAT='($,"Enter the threshold level")' READ,thresh PRINT," " i=FIX(i-1) chgthresh=1 END 'V': i=FIX(i+4) 'W': i=FIX((i-6)>(-1)) 'X': i=FIX(i+9) 'Y': i=FIX((i-11)>(-1)) 'G': BEGIN READ,'Enter image number to skip to: ',go i=FIX(((go-1)>(-1))<(n-2)) END 'H': BEGIN ; tag this image for later PostScript plotting ps_plot=[ps_plot,i] PRINT,'Tagging this image for PostScript plotting.' PRINT,' ' i=FIX(i-1) END 'E': BEGIN HELP,im,im2 STOP,'Type ".con" to resume.' i=FIX(i-1) END 'M': BEGIN ; Display an image slice along a selected meridian IF N_ELEMENTS(merunit) EQ 0 THEN BEGIN ; Open a file to hold information about meridian slices, ; if one is not already open. mername='/jasper/cnsr3_data1/tmp/merdata.tmp' OPENW,merunit,mername,/GET_LUN plotsd='' READ,'Plot s.d. of meridian slice(s)? (Y/[N]) ',plotsd plotsd=(STRUPCASE(plotsd) EQ 'Y') wideslice=1. READ,'Enter meridian slice width (deg; default=1): ',wideslice READ,'Enter the full-scale brightness for the meridian plot: ' $ ,mpfsb ENDIF FOR kmer=0,N_ELEMENTS(double) DO BEGIN ; Only need to select azimuth for 1st image IF kmer EQ 0 THEN BEGIN ; determine image size sim=SIZE(im) simx=sim(1) simy=sim(2) ; position cursor to guide user TVCRS,xsh+simx/2,ysh+61+simy/2,/DEVICE PRINT,'Click on a pixel defining the ' $ +'meridian you want to sample.' ; read cursor position when mouse button pressed CURSOR,xm,ym,/DOWN,/DEVICE TV,REPLICATE(!D.N_COLORS-1,3,3),xm-1,ym-1 ; convert cursor coordinates to image coordinates xm=xm-xsh ym=ym-(ysh+61) pixok=((0 LE xm) AND (xm LE simx) AND $ (0 LE ym) AND (ym LE simy)) IF NOT pixok THEN BEGIN ; Prompt user to enter azimuth from keyboard PRINT,'Pixel outside image!' sunazel,year,month,day,hour,minute,second,sunaz PRINT,'Solar azimuth at image time: ',sunaz,' deg' $ ,FORMAT='(A,F5.1,A)' READ,'Enter the azimuth to sample: ',maz xm=-1 ym=-1 pixok=1B ; Prompt user to return cursor to window so ; MOVCROSS doesn't fail at TVRD. PRINT,'Click in plot window to continue.' CURSOR,jx,jy,/DOWN,/DEVICE ENDIF IF (xm GE 0) AND (ym GE 0) THEN BEGIN PRINT,xm,ym,FORMAT='("xm=",I0," ym=",I0)' ; determine the meridian maz=az(xm,ym) ENDIF PRINT,maz,FORMAT='("Azimuth ",F5.1," deg selected.")' ; Identify the meridian on the image IF KEYWORD_SET(double) $ THEN azval=(nsteps/2-1) $ ELSE azval=!D.N_COLORS-1 waz=WHERE(ABS(az-maz) LT 0.5,nwaz) IF nwaz GT 0 THEN tvim(waz)=azval oaz=(maz+180) MOD 360 waz=WHERE(ABS(oaz-az) LT 0.5,nwaz) IF nwaz GT 0 THEN tvim(waz)=azval TV,tvim,xsh,ysh+61 ; Now save the currently displayed images to a pixmap ; so they can be refreshed later IF N_ELEMENTS(pwin) EQ 0 THEN BEGIN dxsize=!D.X_SIZE dysize=!D.Y_SIZE dwin=!D.WINDOW WINDOW,/FREE,/PIXMAP,XSIZE=dxsize-xsh,YSIZE=dysize pwin=!D.WINDOW WSET,dwin ENDIF WSET,pwin DEVICE,COPY=[xsh,0,dxsize-xsh,dysize,0,0,dwin] WSET,dwin ; Now use routine MERSLICE to return the meridian intensities merslice,ROTATE(im,2),az,za,maz,elev,slice,sd $ ,swidth=wideslice ; Now plot the meridian slice(s) xp=!X.MARGIN*!D.X_CH_SIZE yp=!Y.MARGIN*!D.Y_CH_SIZE posn=[xp(0),yp(0),0.95*xsh-xp(1),!D.Y_SIZE-yp(1)] ;PRINT,'xp,yp,posn' ;PRINT,xp,yp,posn ;PRINT,'!X.S,!Y.S' ;PRINT,!X.S,!Y.S PLOT,[0,1],[0,1] $ ,XTITLE='Elevation Angle (deg)' $ ,XRANGE=[20,160],/XSTYLE $ ,YTITLE='Intensity along Meridian (R)' $ ,YRANGE=[0,mpfsb] $ ,TITLE=STRING(maz,FORMAT='(F5.1)')+' deg azimuth' $ ,/NODATA,POSITION=posn,/DEVICE ; Now refresh the displayed images DEVICE,COPY=[0,0,dxsize-xsh,dysize,xsh,0,pwin] ; Now display the meridian slice pcol=0.6*(!D.N_COLORS-1)/(1.+N_ELEMENTS(double)) OPLOT,elev,slice,COLOR=pcol IF plotsd THEN BEGIN OPLOT,elev,slice+sd,LINESTYLE=1,COLOR=!D.N_COLORS-1 OPLOT,elev,slice-sd,LINESTYLE=1,COLOR=!D.N_COLORS-1 ENDIF IF (xm GE 0) AND (ym GE 0) THEN BEGIN ; Prompt user to return cursor to window so ; MOVCROSS doesn't fail at TVRD. PRINT,'Click in plot window to continue.' CURSOR,jx,jy,/DOWN,/DEVICE ENDIF dmovcross,pkel,pkint,color=pcol ; Write meridian data into the temporary file. WRITEU,merunit,hdb,maz,pkel,pkint $ ,N_ELEMENTS(elev),elev,slice ENDIF ELSE BEGIN ; If a second image is displayed, change the parameters ; and repeat the process for that image. ; First identify the meridian on the second image azval=!D.N_COLORS-1 waz=WHERE(ABS(az2-maz) LT 0.5,nwaz) IF nwaz GT 0 THEN tvim2(waz)=azval oaz=(maz+180) MOD 360 waz=WHERE(ABS(oaz-az2) LT 0.5,nwaz) IF nwaz GT 0 THEN tvim2(waz)=azval TV,tvim2,xsh,ysh2+61 ; Then generate and plot the meridian slice merslice,ROTATE(im2,2),az2,za2,maz,elev,slice,sd $ ,swidth=wideslice pcol=((1.0+0.6)/2.)*(!D.N_COLORS-1) OPLOT,elev,slice,COLOR=pcol dmovcross,pkel,pkint,color=pcol ; Write meridian data into the temporary file. WRITEU,merunit,hdb2,maz,pkel,pkint $ ,N_ELEMENTS(elev),elev,slice ENDELSE ENDFOR READ,'Press to go on.',decis i=FIX(i-1) END 'Q': GOTO,DONE 'Z': BEGIN ; choose a 5 by 5 pixel region to calculate the intensity over rdim(i)=1B maxim=MAX(im,MIN=minim) minim=MAX([0,minim]) tvim=ROTATE(im,2) ; create box around old region - 1 deg about magnetic zenith mzim=INTARR(256,256) mzim(min_c:max_c,min_r)=maxim mzim(min_c:max_c,max_r)=maxim mzim(min_c,min_r:max_r)=maxim mzim(max_c,min_r:max_r)=maxim cursor=1 ; get coordinates of new region (user can choose coordinates from the ; central 64 by 64 pixels) medc=(min_c+max_c)/2 medr=(min_r+max_r)/2 PRINT,medc,medr $ ,FORMAT='("Present region center: (",I0,",",I0,")")' tvsi,(tvim+mzim)<(maxim+1),medc-32,medr-32,medc+31,medr+31 $ ,cursor=cursor,ctable=-1,/noclear ; store the new coordinates ccol(i)=cursor(0) crow(i)=cursor(1) min_c=ccol(i)-2 & max_c=ccol(i)+2 min_r=crow(i)-2 & max_r=crow(i)+2 PRINT,ccol(i),crow(i) $ ,FORMAT='("New region center: (",I0,",",I0,") ",$)' PRINT,az(ccol(i),crow(i)),za(ccol(i),crow(i)) $ ,FORMAT='("(azimuth ",F5.1,", zenith angle ",F4.1,")")' ; calculate the average intensity and the standard deviation in the ; intensity over a 5 by 5 pixel area, centered on the new coordinate zint=STDEV(tvim(min_c:max_c,min_r:max_r),mean) var(i)=zint u(i)=mean mean=STRING(mean,FORMAT='(F5.1)') zint=STRING(zint,FORMAT='(F5.1)') ; put the 5 by 5 pixel box on the primary image tvim(min_c-1:max_c+1,min_r-1)=maxim tvim(min_c-1:max_c+1,max_r+1)=maxim tvim(min_c-1,min_r-1:max_r+1)=maxim tvim(max_c+1,min_r-1:max_r+1)=maxim IF KEYWORD_SET(double) THEN BEGIN bar=BYTSCL(REBIN(BINDGEN(1,256),32,256),MAX=255,MIN=0,TOP=(nsteps/2-1)) tvim=BYTSCL(tvim,MAX=maxim,MIN=minim,TOP=(nsteps/2-1)) ENDIF ELSE BEGIN tvim=BYTSCL(tvim,MAX=maxim,MIN=minim,TOP=(!D.N_COLORS-1)) bar=BYTSCL(REBIN(BINDGEN(1,256),32,256),MAX=255,MIN=0,TOP=(!D.N_COLORS-1)) ENDELSE ; find the zenith and azimuth angles of the new coordinate zen(i)=za(ccol(i),crow(i)) azm(i)=az(ccol(i),crow(i)) ; convert the zenith and azimuth angles to latitude and east longitude angles ; in radians lam1=(90-zen(i))*!DTOR phi1=azm(i)*!DTOR ; calculate the deviation, in degrees, of the new coordinate from the magnetic ; zenith del(i)=!RADEG*ACOS(SIN(lam1)*SIN(lam2)+COS(phi2-phi1)*COS(lam1)*COS(lam2)) szenith=STRING(zen(i),FORMAT='(F4.1)') sazimuth=STRING(azm(i),FORMAT='(F5.1)') sdel=STRING(del(i),FORMAT='(F4.1)') ; display the primary image and results ERASE XYOUTS,xsh,ysh+40 $ ,time+' '+date+' Cam: '+scn+' Filt: '+sfn $ ,/DEVICE,CHARSIZE=csize XYOUTS,xsh,ysh+20 $ ,'Zenith: '+szenith+' Azimuth: '+sazimuth+' Deviation: '+sdel $ ,/DEVICE,CHARSIZE=csize TV,tvim,xsh,ysh+61 PLOT,[0,1],[minim,maxim],/DEVICE,/NODATA $ ,POSITION=[xsh+306,ysh+61,xsh+337,ysh+316] $ ,XSTYLE=4,YRANGE=[minim,maxim],YSTYLE=9,/NOERASE TV,bar,xsh+306,ysh+61 XYOUTS,xsh,ysh $ ,'Zenith intensity: '+mean+' Std.: '+zint $ ,/DEVICE,CHARSIZE=csize IF KEYWORD_SET(double) THEN BEGIN rdim2(i)=1B maxim2=MAX(im2,MIN=minim2) minim2=MAX([0,minim2]) tvim2=ROTATE(im2,2) ; find the zenith angle and azimuth angle that the new pixel chosen on the ; primary image corresponds to on the secondary image. zen2(i)=zen(i) azm2(i)=azm(i) delarr=!RADEG*ACOS(SIN(lam1)*alpha $ +COS((az2*!DTOR)-phi1)*COS(lam1)*beta) cpix=MIN(WHERE(delarr EQ MIN(delarr))) crow2(i)=cpix/256 ; row number on secondary image ccol2(i)=cpix MOD 256 ; column number on secondary image zen2(i)=za2(ccol2(i),crow2(i)) ; zenith angle of new pixel azm2(i)=az2(ccol2(i),crow2(i)) ; azimuth angle of new pixel lam12=(90-zen2(i))*!DTOR phi12=azm2(i)*!DTOR min_c2=ccol2(i)-2 & max_c2=ccol2(i)+2 min_r2=crow2(i)-2 & max_r2=crow2(i)+2 ; calculate the intensity over the 5 by 5 pixel box centered on the ; 'new pixel' zint=STDEV(tvim2(min_c2:max_c2,min_r2:max_r2),mean) var2(i)=zint u2(i)=mean mean=STRING(mean,FORMAT='(F5.1)') zint=STRING(zint,FORMAT='(F5.1)') ; put the 5 by 5 pixel box on the secondary image tvim2(min_c2-1:max_c2+1,min_r2-1)=maxim2 tvim2(min_c2-1:max_c2+1,max_r2+1)=maxim2 tvim2(min_c2-1,min_r2-1:max_r2+1)=maxim2 tvim2(max_c2+1,min_r2-1:max_r2+1)=maxim2 bar2=BYTSCL(REBIN(BINDGEN(1,256),32,256) $ ,MAX=255,MIN=0,TOP=(nsteps/2-1)) tvim2=BYTSCL(tvim2,MAX=maxim,MIN=minim,TOP=(nsteps/2-1)) bar2=bar2+BYTE(nsteps/2+1) tvim2=tvim2+BYTE(nsteps/2+1) ; calculate the deviation, in degrees, of the new coordinate from the magnetic ; zenith (for the secondary image) del2(i)=!RADEG*ACOS(SIN(lam12)*SIN(lam22) $ +COS(phi22-phi12)*COS(lam12)*COS(lam22)) szenith=STRING(zen2(i),FORMAT='(F4.1)') sazimuth=STRING(azm2(i),FORMAT='(F5.1)') sdel=STRING(del2(i),FORMAT='(F4.1)') ; display the secondary image and results XYOUTS,xsh,ysh2+40 $ ,time+' '+date+' Cam2: '+scn2+' Filt: '+sfn2 $ ,/DEVICE,CHARSIZE=csize XYOUTS,xsh,ysh2+20 $ ,'Zenith: '+szenith+' Azimuth: '+sazimuth+' Deviation: '+sdel $ ,/DEVICE,CHARSIZE=csize TV,tvim2,xsh,ysh2+61 PLOT,[0,1],[minim2,maxim2],/DEVICE,/NODATA $ ,POSITION=[xsh+306,ysh2+61,xsh+337,ysh2+316] $ ,XSTYLE=4,YRANGE=[minim2,maxim2],YSTYLE=9,/NOERASE TV,bar2,xsh+306,ysh2+61 XYOUTS,xsh,ysh2 $ ,'Zenith intensity: '+mean+' Std.: '+zint $ ,/DEVICE,CHARSIZE=csize ENDIF i=FIX(i-1) norefresh=1B END ELSE: ENDCASE ENDREP UNTIL 0 EQ 1 DONE: BEGIN PRINT,'Done; tying up loose ends.' IF N_ELEMENTS(merunit) GT 0 THEN BEGIN ; Close the temporary file holding meridian data. FREE_LUN,merunit ; Plot the selected meridian slices to a PostScript file. psmslice,year,month,day,mpfsb ENDIF IF N_ELEMENTS(ps_plot) GT 1 THEN BEGIN ; There are images to plot out, too. ps_plot=ps_plot(1:*) ps_plot=ps_plot(uniq(ps_plot(SORT(ps_plot)))) PRINT,'Outputting the following images to a PostScript file:' PRINT,STRTRIM(ps_plot,2) bwc='' REPEAT BEGIN READ,'Do you want a black-and-white or colour plot? [B/C] ',bwc bwc=STRUPCASE(bwc) ENDREP UNTIL (bwc EQ 'B') OR (bwc EQ 'C') psname='' READ,'Enter a name for the PostScript file: ',psname ; Open PostScript file for the type of output requested by the user. CASE bwc OF 'B': BEGIN psopen,psname DEVICE,XSIZE=8.143,YSIZE=10.0,XOFFSET=0.1785,YOFFSET=1.75 $ ,/INCHES END 'C': BEGIN psopen,psname,/COLOR,BITS_PER_PIXEL=8 DEVICE,XSIZE=8.143,YSIZE=10.0,XOFFSET=0.1785,YOFFSET=1.75 $ ,/INCHES END ENDCASE kmax=N_ELEMENTS(ps_plot) IF KEYWORD_SET(double) THEN kmax=2*kmax ; set up the images so that there will be 9 images per page - 3 images ; per row and 3 rows of images. FOR k=0,kmax-1 DO BEGIN s=(k MOD 9)/3 e=k MOD 3 ; if two images were displayed, then images will be saved alternating ; between primary and secondary images across the page IF KEYWORD_SET(double) THEN BEGIN i=ps_plot(k/2) IF (k MOD 2) NE 0 THEN BEGIN j=keep2(i) im=BYTE(files2(j),512,256,256) hdb=BYTE(files2(j),0,512) dev=del2(i) std=var2(i) mean=u2(i) azimuth=azm2(i) zenith=zen2(i) col=ccol2(i) row=crow2(i) region=rdim2(i) scn=cn2 sfn=fn2 mic=minc2 mac=maxc2 mir=minr2 mar=maxr2 ENDIF ELSE BEGIN j=keep(i) im=BYTE(files(j),512,256,256) hdb=BYTE(files(j),0,512) dev=del(i) std=var(i) mean=u(i) azimuth=azm(i) zenith=zen(i) col=ccol(i) row=crow(i) region=rdim(i) scn=cn sfn=fn mic=minc mac=maxc mir=minr mar=maxr ENDELSE j=keep(i) hour=utsec(j)/3600L minute=(utsec(j)-hour*3600L)/60L second=utsec(j)-hour*3600L-minute*60L time=STRING(hour,minute,second $ ,FORMAT='(I2.2,":",I2.2,":",I2.2)') hd=gethd(hdb) power=hd.proc.base2 scale=hd.proc.scale_factor IF scale GT 0 THEN BEGIN IF power EQ 1 THEN im=LONG(im)*(2L^scale) IF power EQ 0 THEN im=LONG(im)*scale ENDIF IF level(i) GT 0 THEN im=im