;Program to read in report files generated by PAHFITS and produce ;Countour plots pro maptoom17 COMMON FILEBLOCK, allreportfiles,reportcount,setplot,source,image set_plot,'x' setplot='X' ; CK -- These device settings are for when plotting on the terminal window. device, true_color=24 device, decompose=0,retain=2 ;print, 'WARNING: currently set to plot to postscript file! ' ;set_plot, 'ps' ;setplot='ps' ;device, file='h2_s2map.eps', bits_per_pixel=8, /color, /encapsulated,$ ; /landscape ;Top level path to report files ;path = '/peto6/claudia/PDRs/TEST/' ;path = '/Users/mark/spitzer/m17sw' path = '/home/mwolfire/spitzer/m17sw' image = "/home/mwolfire/spitzer/tools/remy/M17-UC1.I4.rd.fits" ;Select Directory for report files reportpath = DIALOG_PICKFILE(path=path,/DIRECTORY,title='Select report path') ;Search for ALL report files generatred by PAHFITS allreportfiles = FILE_SEARCH(reportpath,'report*',COUNT=reportcount) print,reportcount,' report files found' ;Start description for CW_FORM desc_disp = [ $ '0, LABEL, Select Lines to Map, CENTER', $ '1,BASE,,COLUMN,FRAME', $ '0, BUTTON, [S IV] |H_2S(2) |[Ne II],ROW,TAG=LINESELECT1', $ '0, BUTTON, [Ne III]|H_2S(1) |[S III],ROW,TAG=LINESELECT2', $ '0, BUTTON, [O IV] |[Fe II] |H_2S(0),ROW,TAG=LINESELECT3', $ '2, BUTTON, [S III] | [Si II]| ,ROW,TAG=LINESELECT4', $ '1,BASE,,ROW,FRAME',$ '0,BUTTON, Ratio: |Contour:,ROW,TAG=RATCON, EXCLUSIVE, SET_VALUE=[0]',$ '2, DROPLIST, None |[S IV] |H_2S(2) |[Ne II] |[Ne III] |H_2S(1) |[S III]'+$ '|[O IV]|[Fe II] |H_2S(0) |[S III] |[Si II],'+$ 'ROW,TAG=LINESELECTr1', $ '1,BASE,,ROW,FRAME',$ '2,BUTTON, Yes: |No:,ROW,TAG=DSTRIP, EXCLUSIVE, SET_VALUE=[0],'+$ 'LABEL_LEFT=DeStripe:',$ '1,BASE,,ROW,FRAME',$ '0,BUTTON,Yes:|No:,ROW,TAG=FULLR, EXCLUSIVE, SET_VALUE=[0],'+$ 'LABEL_LEFT=Full Rangee:',$ '2,FLOAT, 0, TAG=MAXPLOTV',$ '1,BASE,,ROW,FRAME',$ '0,BUTTON,Mapit!, TAG=MAPIT', $ '0,BUTTON,Close Maps, TAG=CLOSEMAPS', $ '0,BUTTON,Close All, TAG=CLOSEALL', $ '2,BUTTON, Quit, Quit, TAG=Quit'] base_disp = WIDGET_BASE(TITLE='Maptoo Control Panel',/COLUMN) b_disp = CW_FORM(base_disp,desc_disp,/COLUMN,UNAME='maps') ;Rchisaq slider is not a FORMS widget b_slider = CW_FSLIDER(base_disp,MAXIMUM = 100, SCROLL=1, $ MINIMUM=1,/DRAG,VALUE=50,UNAME='rchisq', $ TITLE='Select Maximum Reduced Chisq',FORMAT='(5X ,I3)') WIDGET_CONTROL,base_disp,/REAL XMANAGER,'maptoo',base_disp,/NO_BLOCK END ;End Widget ------------------------------------------------------ ;Start Event Handler --------------------------------------------- pro maptoo_event,event COMMON FILEBLOCK, allreportfiles,reportcount,setplot,source,image WIDGET_CONTROL,event.id, GET_VALUE=value uname=WIDGET_INFO(event.id,/UNAME) ;rchisq is not a forms widget, do nothing when sliding IF uname EQ 'rchisq' THEN RETURN ;Forms Widget IF uname EQ 'maps' THEN BEGIN tag = event.tag ; print, tag IF event.quit EQ 1 THEN BEGIN WIDGET_CONTROL,event.top,/DESTROY ENDIF IF tag EQ 'LINESELECT1' OR tag EQ 'LINESELECT2' OR $ tag EQ 'LINESELECT3' OR tag EQ 'LINESELECT4' THEN RETURN IF tag EQ 'CLOSEALL' THEN BEGIN FOR loop = 1, 50 DO BEGIN IF !D.WINDOW NE -1 THEN WDELETE, !D.WINDOW ENDFOR ENDIF ;How to close maps but leave up chisq? tag = cloasmaps ; IF tag EQ 'CLOSEMAPS' THEN BEGIN ; FOR loop = 1, 20 DO BEGIN ; IF !D.WINDOW NE -1 AND !D.WINDOW NE 0 THEN WDELETE, !D.WINDOW ; ENDFOR ; ENDIF IF tag EQ 'MAPIT' THEN BEGIN print,value.lineselect1 print,value.lineselect2 print,value.lineselect3 print,value.lineselect4 ; print, tag nmaps = TOTAL(value.lineselect1) + $ TOTAL(value.lineselect2) + $ TOTAL(value.lineselect3) + $ TOTAL(value.lineselect4) linearray = [value.lineselect1, value.lineselect2, $ value.lineselect3, value.lineselect4] print, 'line array',linearray nsh = TOTAL(value.lineselect1)+TOTAL(value.lineselect2) nlh = TOTAL(value.lineselect3)+TOTAL(value.lineselect4) IF nmaps EQ 0 THEN BEGIN nomap = DIALOG_MESSAGE('Please select a line to map') RETURN ENDIF widrchisq = WIDGET_INFO(event.top,FIND_BY_UNAME='rchisq') WIDGET_CONTROL,widrchisq,GET_VALUE=rchisqlim print,rchisqlim ;Search strings in report files mapsearchstring = ['*SIV*','*H2 S(2)*','*NeII*', $ '*NeIII*','*H2 S(1)*','*SIII*18 lam*', $ '*OIV*','*FeII*','*H2 S(0)*', $ '*SIII*33 lam*','*SiII*'] sivs = textoidl('[S IV] 10.51 \mum') h2s2s = textoidl('H_2 S(2) 12.28 \mum') neiis = textoidl('[Ne II] 12.81 \mum') neiiis = textoidl('[Ne III] 15.55 \mum') h2s1s = textoidl('H_2 S(1) 17.03 \mum') siiis = textoidl('[S III] 18.71 \mum') oivs = textoidl('[O IV] 25.91 \mum') feiis = textoidl('[Fe II] 25.99 \mum') h2s0s = textoidl('H_2 S(0) 28.22 \mum') siiis2 = textoidl('[S III] 33.48 \mum') si2s = textoidl('[Si II] 34.81 \mum') maptitlestring = [sivs, h2s2s,neiis,$ neiiis, h2s1s,siiis,$ oivs, feiis,h2s0s,$ siiis2,si2s] lineindex = WHERE(linearray EQ 1) shorthigh = 0 longhigh = 0 ;Read in reports and get position and line intensity FOR nfiles = 0,reportcount-1 DO BEGIN reportfile = STRARR(200) reportfiles = ' ' ;Get y position of report to skip rows for de stripe ypos = STRPOS(allreportfiles[nfiles],'y') ylen=STRLEN(allreportfiles[nfiles])-ypos-1 ypx = float(STRMID(allreportfiles[nfiles],ypos+1,ylen)) ;AOR test to join AORs between SH-06 and SH-07 and ; LH-02 and LH-03 aorpos = STRPOS(allreportfiles[nfiles],'-PDR-',/REVERSE_SEARCH) aortest = (STRMID(allreportfiles[nfiles],aorpos+5,5)) modtest = (STRMID(aortest,0,2)) m17shy = [0,2,4,6] shall = [0,1,2,3,4,6,7] m17shjoin = [0,2,4,6,7] m17lhy = [0,1,2,4,5,7,8,9,10,12,13] lhall = [0,1,2,4,5,7,8,9,10,12,13] m17lhjoin = [0,1,2,4,5,7,8,9,10,12,13] ;Skip every other line (needs fixing for LH) IF modtest EQ 'SH' THEN BEGIN IF value.dstrip EQ 0 AND aortest NE 'SH-06' THEN BEGIN ystrip = m17shy ENDIF ELSE BEGIN ystrip = m17shjoin ENDELSE ;Do all lines IF value.dstrip EQ 1 THEN ystrip = shall ENDIF IF modtest EQ 'LH' THEN BEGIN IF value.dstrip EQ 0 AND aortest NE 'LH-03' THEN BEGIN ystrip = m17lhy ENDIF ELSE BEGIN ystrip = m17lhjoin ENDELSE ;Do all lines IF value.dstrip EQ 1 THEN ystrip = lhall ENDIF IF WHERE(ypx EQ ystrip) NE -1 THEN BEGIN OPENR,unit,allreportfiles[nfiles],/GET_LUN i = 0 WHILE NOT EOF(unit) DO BEGIN READF,unit,reportfiles ;Load each line of reportfile into an array reportfile[i] = reportfiles IF i LT (200-1) THEN i = i+1 ELSE STOP, 'Too many lines' ENDWHILE CLOSE,unit FREE_LUN,unit ;Truncate file reportfile = reportfile[0:i-1] ;Extract AOR from string aorpos = STRPOS(allreportfiles[nfiles],'-PDR-',/REVERSE_SEARCH) aor = (STRMID(allreportfiles[nfiles],aorpos+5,5)) ;Alpha and dec are stored on last line of report file ; but are differerent for each module. ;Search for Line intensity; check that SH and/or LH data are present ;Extract all line fluxes at once IF STRMATCH(allreportfiles[nfiles],'*/ch1/*') THEN BEGIN shorthigh = 1 IF N_ELEMENTS(fluxsiv) EQ 0 THEN BEGIN fluxsiv =(STRSPLIT(reportfile[45],' ',/EXTRACT))[6] fluxh2s2 =(STRSPLIT(reportfile[48],' ',/EXTRACT))[6] fluxneiii =(STRSPLIT(reportfile[51],' ',/EXTRACT))[6] fluxneii=(STRSPLIT(reportfile[54],' ',/EXTRACT))[6] fluxh2s1 =(STRSPLIT(reportfile[57],' ',/EXTRACT))[6] fluxsiii =(STRSPLIT(reportfile[60],' ',/EXTRACT))[6] alphash = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[2] decsh = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[3] aorsh = aor rchisqsh = (STRSPLIT(reportfile[4],' ',/EXTRACT))[5] ;Check for *** in rchisqsh and set to 500 tmatch = STRMATCH(rchisqsh,"*\*\*\*\**") if tmatch EQ '1' THEN BEGIN rchisqsh = 500 PRINT,'rchisqsh not defined: ',$ allreportfiles[nfiles] ENDIF ENDIF ELSE BEGIN fluxsiv0 =(STRSPLIT(reportfile[45],' ',/EXTRACT))[6] fluxh2s20 =(STRSPLIT(reportfile[48],' ',/EXTRACT))[6] fluxneiii0 =(STRSPLIT(reportfile[51],' ',/EXTRACT))[6] fluxneii0=(STRSPLIT(reportfile[54],' ',/EXTRACT))[6] fluxh2s10 =(STRSPLIT(reportfile[57],' ',/EXTRACT))[6] fluxsiii0 =(STRSPLIT(reportfile[60],' ',/EXTRACT))[6] alphash0 = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[2] decsh0 = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[3] aorsh0 = aor rchisqsh0 = (STRSPLIT(reportfile[4],' ',/EXTRACT))[5] ;Check for *** in rchisqsh0 and set to 500 tmatch = STRMATCH(rchisqsh0,"*\*\*\*\**") if tmatch EQ '1' THEN BEGIN rchisqsh0 = 500 PRINT,'rchisqsh0 not defined: ',$ allreportfiles[nfiles] ENDIF ENDELSE ENDIF IF STRMATCH(allreportfiles[nfiles],'*/ch3/*') THEN BEGIN longhigh = 1 IF N_ELEMENTS(fluxoiv) EQ 0 THEN BEGIN fluxoiv =(STRSPLIT(reportfile[63],' ',/EXTRACT))[6] fluxfeii =(STRSPLIT(reportfile[66],' ',/EXTRACT))[6] fluxh2s0 =(STRSPLIT(reportfile[69],' ',/EXTRACT))[6] fluxsiii2 =(STRSPLIT(reportfile[72],' ',/EXTRACT))[6] fluxsi2 =(STRSPLIT(reportfile[75],' ',/EXTRACT))[6] alphalh = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[2] declh = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[3] aorlh = aor rchisqlh = (STRSPLIT(reportfile[4],' ',/EXTRACT))[5] ;Check for *** in rchisqlh and set to 500 tmatch = STRMATCH(rchisqlh,"*\*\*\*\**") if tmatch EQ '1' THEN BEGIN rchisqlh = 500 PRINT,'rchisqlh not defined: ',$ allreportfiles[nfiles] ENDIF ENDIF ELSE BEGIN fluxoiv0 =(STRSPLIT(reportfile[63],' ',/EXTRACT))[6] fluxfeii0 =(STRSPLIT(reportfile[66],' ',/EXTRACT))[6] fluxh2s00 =(STRSPLIT(reportfile[69],' ',/EXTRACT))[6] fluxsiii20 =(STRSPLIT(reportfile[72],' ',/EXTRACT))[6] fluxsi20 =(STRSPLIT(reportfile[75],' ',/EXTRACT))[6] alphalh0 = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[2] declh0 = (STRSPLIT(reportfile[i-1],' ',/EXTRACT))[3] aorlh0 = aor rchisqlh0 = (STRSPLIT(reportfile[4],' ',/EXTRACT))[5] ;Check for *** in rchisqlh0 and set to 500 tmatch = STRMATCH(rchisqlh0,"*\*\*\*\**") if tmatch EQ '1' THEN BEGIN rchisqlh0 = 500 PRINT,'rchisqlh0 not defined: ',$ allreportfiles[nfiles] ENDIF ENDELSE ENDIF IF nfiles GT 0 THEN BEGIN ;Check for SH report files IF STRMATCH(allreportfiles[nfiles],'*/ch1/*') AND $ N_ELEMENTS(fluxsiv0) GT 0 THEN BEGIN fluxsiv = [fluxsiv,fluxsiv0] fluxh2s2 = [fluxh2s2,fluxh2s20] fluxneiii = [fluxneiii,fluxneiii0] fluxneii = [fluxneii,fluxneii0] fluxh2s1 = [fluxh2s1,fluxh2s10] fluxsiii = [fluxsiii,fluxsiii0] alphash = [alphash,alphash0] decsh = [decsh,decsh0] aorsh = [aorsh,aorsh0] rchisqsh =[rchisqsh,rchisqsh0] ENDIF ;Check for LH report files IF STRMATCH(allreportfiles[nfiles],'*/ch3/*') AND $ N_ELEMENTS(fluxoiv0) GT 0 THEN BEGIN fluxoiv = [fluxoiv,fluxoiv0] fluxfeii = [fluxfeii,fluxfeii0] fluxh2s0= [fluxh2s0,fluxh2s00] fluxsiii2 = [fluxsiii2,fluxsiii20] fluxsi2 = [fluxsi2,fluxsi20] alphalh = [alphalh,alphalh0] declh = [declh,declh0] aorlh = [aorlh,aorlh0] rchisqlh =[rchisqlh,rchisqlh0] ENDIF ENDIF ENDIF ENDFOR ;end reading in report files loadct,39,/silent, ncolors=255 IF nsh GT 0 AND shorthigh EQ 0 THEN BEGIN oops = DIALOG_MESSAGE(' No SH Data ') RETURN ENDIF IF shorthigh EQ 1 THEN BEGIN IF setplot EQ 'X' THEN BEGIN WINDOW,/free WSET,!D.WINDOW ENDIF ;Plot maps of reduced Chi^2 SH xtitle=textoidl('SH Reduced \chi^2') rchmaxplt = 600.0 nlev = 255 ulev=(findgen(nlev+1)/nlev) *$ ((rchmaxplt)-min(float(rchisqsh)))+min(float(rchisqsh)) ; plot,HISTOGRAM(float(rchisqsh)),XTITLE=xtitle,$ ; YRANGE=[0,20],PSYM=10 contour,float(rchisqsh),float(alphash),float(decsh),levels=ulev,$ XRANGE=[max(float(alphash)),min(float(alphash))],$ yr=[min(float(decsh)),max(float(decsh))],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=textoidl('SH Reduced \chi^2'), $ POS=[0.15,0.15,0.75,0.90], $ c_colors = Indgen(nlev), charsize=1.1 ; colorbar,RANGE=[min(float(rchisqsh)),max(float(rchisqsh))],$ colorbar,RANGE=[min(float(rchisqsh)),rchmaxplt],$ /VERTICAL,FORMAT='(E8.1)',/LEFT, $ ncolors=255, charsize=1.1 ENDIF IF nlh GT 0 AND longhigh EQ 0 THEN BEGIN oops = DIALOG_MESSAGE(' No LH Data ') RETURN ENDIF IF longhigh EQ 1 THEN BEGIN IF setplot EQ 'X' THEN BEGIN WINDOW,/free WSET,!D.WINDOW ENDIF ;Plot maps of reduced Chi^2 LH rchmaxplt = 600.0 nlev = 255 ulev=(findgen(nlev+1)/nlev) *$ ((rchmaxplt)-min(float(rchisqlh)))+min(float(rchisqlh)) contour,float(rchisqlh),float(alphalh),float(declh),levels=ulev,$ XRANGE=[max(float(alphalh)),min(float(alphalh))],$ yr=[min(float(declh)),max(float(declh))],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=textoidl('LH Reduced \chi^2'), $ POS=[0.15,0.15,0.75,0.90], $ c_colors = Indgen(nlev), charsize=1.1 colorbar,RANGE=[min(float(rchisqlh)),rchmaxplt],$ /VERTICAL,FORMAT='(E8.1)',/LEFT, $ ncolors=255, charsize=1.1 ENDIF FOR imap = 0,nmaps-1 DO BEGIN linestring = mapsearchstring[lineindex[imap]] maptitle = maptitlestring[lineindex[imap]] CASE lineindex[imap] of 0: BEGIN flux = float(fluxsiv) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 1: BEGIN flux = float(fluxh2s2) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 2: BEGIN flux = float(fluxneii) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 3: BEGIN flux = float(fluxneiii) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 4: BEGIN flux = float(fluxh2s1) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 5: BEGIN flux = float(fluxsiii) alpha = alphash dec = decsh aorplt = aorsh rchisq = float(rchisqsh) END 6: BEGIN flux = float(fluxoiv) alpha = alphalh dec = declh aorplt = aorlh rchisq = float(rchisqlh) END 7: BEGIN flux = float(fluxfeii) alpha = alphalh dec = declh aorplt = aorlh rchisq = float(rchisqlh) END 8: BEGIN flux = float(fluxh2s0) alpha = alphalh dec = declh aorplt = aorlh rchisq = float(rchisqlh) END 9: BEGIN flux = float(fluxsiii2) alpha = alphalh dec = declh aorplt = aorlh rchisq = float(rchisqlh) END 10: BEGIN flux = float(fluxsi2) alpha = alphalh dec = declh aorplt = aorlh rchisq = float(rchisqlh) END ENDCASE IF(N_ELEMENTS(flux) EQ 1) THEN BEGIN print, 'no flux' RETURN ENDIF alpha= (float(alpha)) dec = (float(dec)) rchisq = float(rchisq) ;Extract points less than maximum iplotr1 = where(rchisq LT rchisqlim AND $ (aorplt EQ 'SH-01' OR aorplt EQ 'SH-02'$ OR aorplt EQ 'SH-03' OR aorplt EQ 'SH-04'$ OR aorplt EQ 'SH-05' OR aorplt EQ 'SH-06'$ OR aorplt EQ 'LH-01A' OR aorplt EQ 'LH-02'$ )) iplotr2 = where(rchisq LT rchisqlim AND $ (aorplt EQ 'SH-07' OR aorplt EQ 'SH-08'$ OR aorplt EQ 'SH-09' OR aorplt EQ 'SH-10'$ OR aorplt EQ 'SH-11'$ OR aorplt EQ 'LH-03' OR aorplt EQ 'LH-04'$ )) iplot = [iplotr1,iplotr2] maxf = max(flux[iplot]) minf = min(flux[iplot]) ; Use input maxf value IF value.fullr EQ 1 AND value.maxplotv NE 0 THEN BEGIN maxf = value.maxplotv minf = 0.0 ENDIF IF imap EQ 0 THEN BEGIN maxfden = maxf minfden = minf ENDIF ;maxf = 1.7e-6 S(2) ;maxf = 4.0e-7 S(1) ;maxf = 3.0e-7 S(0) ;hist = HIST_EQUAL(flux[iplot]) ;histsum = TOTAL(hist) ;histscale = hist*(histsum/nlev) ;print,hist[UNIQ(hist[(SORT((hist)))])] ;ulev = hist[UNIQ(hist[SORT((hist))])]*(maxf-minf)+minf maxra = max(alpha[iplot]) minra = min(alpha[iplot]) maxdec= max(dec[iplot]) mindec= min(dec[iplot]) !x.style=1 ; CK -- setting these to 1 forces the xrange to be ; the specified range !y.style=1 ; same goes for yrange IF setplot EQ 'X' THEN BEGIN WINDOW,/free WSET,!D.WINDOW ENDIF ;PLOT Box ;print,'iplot ',iplot ;print,'iplotr1 ',iplotr1 ;print,'iplotr2 ',iplotr2 ;------------------------------------------------------------- ;Start Contour overlays on fits image ; loadct,0,/silent loadct,39,/silent, ncolors=255 im=mrdfits(image,0,hdr) > 0.1 ; for m17, extract an even smaller image hextract,im,hdr,225,425,125,325 ; extents of image extast,hdr,ast s=size(im) xx=[s[1],0] yy=[0,s[2]] if strmid(sxpar(hdr,'CTYPE1'),0,2) eq 'RA' then begin xy2ad,xx,yy,ast,ra_ra,de_ra end else begin xy2ad,xx,yy,ast,lra,bra euler,ra_ra,de_ra,lra,bra,2 end ra = alpha[iplotr1] de = dec[iplotr1] ; set output grid spacing to be ~1/4 of the original grid ngrid=4*sqrt(1.*n_elements(ra)) bin=(max(ra)-min(ra))/ngrid triangulate,ra,de,tr,b linemap=trigrid(ra,de,flux[iplotr1],tr,[bin,bin],$ [ra_ra[0],de_ra[0],ra_ra[1],de_ra[1]]) linemap=reverse(linemap,1) imcontour,im,hdr,/type,levels=[0] oplotimage,bytscl(im) loadct,39,/silent, ncolors=255 nlev=12 ulev=(findgen(nlev+1)/nlev) *(maxf-minf)+minf sl = size(linemap) nx=sl[1] ny=sl[2] print,'ulev ',ulev ;Plot first set of AORs contour,linemap,$ findgen(nx)*s[1]/(nx-1),findgen(ny)/(ny-1)*s[2], $ levels=ulev,$ XTITLE='!17Right Ascension',YTITLE='Declination',$ TITLE=maptitle, C_COLORS=[200], $ charsize=1.1,/overplot,/closed ;Plot snd set of AORs ra = alpha[iplotr2] de = dec[iplotr2] ; set output grid spacing to be ~1/4 of the original grid ngrid=4*sqrt(1.*n_elements(ra)) bin=(max(ra)-min(ra))/ngrid triangulate,ra,de,tr,b linemap=trigrid(ra,de,flux[iplotr2],tr,[bin,bin],$ [ra_ra[0],de_ra[0],ra_ra[1],de_ra[1]]) linemap=reverse(linemap,1) sl = size(linemap) nx=sl[1] ny=sl[2] contour,linemap, $ findgen(nx)*s[1]/(nx-1),findgen(ny)/(ny-1)*s[2], $ levels=ulev,$ XTITLE='!17Right Ascension',YTITLE='Declination',$ TITLE=maptitle, $ c_colors = [200], charsize=1.1,/overplot,closed=0 ; End contour overlays on fits image ;------------------------------------------------------------------ IF setplot EQ 'X' THEN BEGIN WINDOW,/free WSET,!D.WINDOW ENDIF nlev=255 ; linear levels - change this to suit ulev=(findgen(nlev+1)/nlev) *(maxf-minf)+minf contour,flux[iplot],alpha[iplot],dec[iplot],levels=ulev,$ XRANGE=[maxra,minra],yr=[mindec,maxdec],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitle, $ POS=[0.15,0.15,0.75,0.90], $ c_colors = Indgen(nlev), charsize=1.1,/NODATA ;Plot first set of AORs contour,flux[iplotr1],alpha[iplotr1],dec[iplotr1],levels=ulev,$ XRANGE=[maxra,minra],yr=[mindec,maxdec],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitle, $ POS=[0.15,0.15,0.75,0.90], $ c_colors = Indgen(nlev), charsize=1.1,/NOERASE ;Plot 2nd set of AORs contour,flux[iplotr2],alpha[iplotr2],dec[iplotr2],levels=ulev,$ XRANGE=[maxra,minra],yr=[mindec,maxdec],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitle, $ POS=[0.15,0.15,0.75,0.90], $ c_colors = Indgen(nlev), charsize=1.1,/NOERASE colorbar,RANGE=[minf,maxf],/VERTICAL,FORMAT='(E8.1)',/LEFT, $ ncolors=255, charsize=1.1 print,'ratcon ',value.ratcon print,'lineselectr1 ',value.lineselectr1 IF value.lineselectr1 NE 0 THEN BEGIN ;Ratio or contour map ;Check that LH data are found IF value.lineselectr1 GE 7 AND longhigh EQ 0 THEN BEGIN oops = DIALOG_MESSAGE('Long High Data not Found') RETURN ENDIF ;Check that SH data are found IF value.lineselectr1 LT 7 AND shorthigh EQ 0 THEN BEGIN oops = DIALOG_MESSAGE('Short High Data not Found') RETURN ENDIF CASE value.lineselectr1-1 of 0: BEGIN fluxr = float(fluxsiv) alpha = alphash dec = decsh aorplt = aorsh END 1: BEGIN fluxr = float(fluxh2s2) ; alpha = alphash ; dec = decsh ; aorplt = aorsh END 2: BEGIN fluxr = float(fluxneii) alpha = alphash dec = decsh aorplt = aorsh END 3: BEGIN fluxr = float(fluxneiii) alpha = alphash dec = decsh aorplt = aorsh END 4: BEGIN fluxr = float(fluxh2s1) ; alpha = alphash ; dec = decsh ; aorplt = aorsh END 5: BEGIN fluxr = float(fluxsiii) alpha = alphash dec = decsh aorplt = aorsh END 6: BEGIN fluxr = float(fluxoiv) alpha = alphalh dec = declh aorplt = aorlh END 7: BEGIN fluxr = float(fluxfeii) alpha = alphalh dec = declh aorplt = aorlh END 8: BEGIN fluxr = float(fluxh2s0) alpha = alphalh dec = declh aorplt = aorlh END 9: BEGIN fluxr = float(fluxsiii2) alpha = alphalh dec = declh aorplt = aorlh END 10: BEGIN fluxr = float(fluxsi2) alpha = alphalh dec = declh aorplt = aorlh END ENDCASE IF value.ratcon EQ 0 THEN BEGIN ; Ratio map ;nlev=10 ;ulev=(findgen(nlev+1)/nlev)^2 *(maxfden-minfden)+minfden print, N_ELEMENTS(fluxr) print, 'flux',(flux[iplot]) ; ratioplt = (fluxr[iplot]/flux[iplot]) > 0 ; ratioplt = (fluxr[iplot]/flux[iplot]) < 5.0e1 ratioplt = (fluxr/flux) > 0 ratioplt = (fluxr/flux) < 5.0e1 maxf = max(ratioplt) minf = min(ratioplt) maxf = 5.0 print,'maxf ',maxf, 'minf ',minf IF maxf NE minf THEN BEGIN ulev=(findgen(nlev+1)/nlev) *(maxf-minf)+minf IF setplot EQ 'X' THEN BEGIN WINDOW,/FREE WSET,!D.WINDOW ENDIF maptitler = $ maptitlestring[value.lineselectr1-1]+'/'+maptitle print, ratioplt[iplot] contour,ratioplt[iplot],alpha[iplot],dec[iplot],$ levels=ulev,XRANGE=[maxra,minra],$ yr=[mindec,maxdec], $ /FILL,XTITLE='!17Right Ascension',$ YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitler, $ POS=[0.15,0.15,0.75,0.90],$ c_colors = Indgen(nlev), charsize=1.1,/NODATA contour,ratioplt[iplotr1],alpha[iplotr1],dec[iplotr1],$ levels=ulev,XRANGE=[maxra,minra],$ yr=[mindec,maxdec], $ /FILL,XTITLE='!17Right Ascension',$ YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitler, $ POS=[0.15,0.15,0.75,0.90],$ c_colors = Indgen(nlev), charsize=1.1,/NOERASE contour,ratioplt[iplotr2],alpha[iplotr2],dec[iplotr2],$ levels=ulev,XRANGE=[maxra,minra],$ yr=[mindec,maxdec], $ /FILL,XTITLE='!17Right Ascension',$ YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitler, $ POS=[0.15,0.15,0.75,0.90],$ c_colors = Indgen(nlev), charsize=1.1,/NOERASE colorbar,RANGE=[minf,maxf],$ /VERTICAL,FORMAT='(E8.1)',/LEFT, $ ncolors=255, charsize=1.1 ENDIF ENDIF IF value.ratcon EQ 1 THEN BEGIN ;Contour plot IF setplot EQ 'X' THEN BEGIN WINDOW,/FREE WSET,!D.WINDOW ; window, imap+20 ; wset,imap+20 ENDIF maptitler = $ maptitlestring[value.lineselectr1-1]+' on '+maptitle contour,flux[iplot],alpha[iplot],dec[iplot],$ levels=ulev,XRANGE=[maxra,minra],$ yr=[mindec,maxdec],$ /FILL,XTITLE='!17Right Ascension',YTITLE='Declination',$ XTICKFORMAT='radecformat',YTICKFORMAT='radecformat',$ /IRREGULAR, TITLE=maptitler, $ POS=[0.15,0.15,0.75,0.90],$ c_colors = Indgen(nlev), charsize=1.1 colorbar,RANGE=[minf,maxf],$ /VERTICAL,FORMAT='(E8.1)',/LEFT, $ ncolors=255, charsize=1.1 maxf = max(fluxr[iplot]) minf = min(fluxr[iplot]) nlevr = 10 ulev=(findgen(nlevr+1)/nlevr)^2 *(maxf-minf)+minf contour,fluxr[iplot],alpha[iplot],dec[iplot], $ levels=ulev,XRANGE=[maxra,minra],$ yrange=[mindec,maxdec],$ /IRREGULAR, /OVERPLOT, $ POS=[0.15,0.15,0.75,0.90],$ charsize=1.1 ENDIF ENDIF ;device,/close set_plot, 'x' ENDFOR ENDIF ENDIF END function radecformat,axis,index,value ;Function to return ra, dec in hex ; CK - made changes to manipulate the tick labels. Currently ; set to write RA in hr min sec for 1st label and only sec after that. ; For Dec, the top label is set to print, deg min as and other labels only ; min as. However, this will only work if there are 7 labels as I had ; to hardwire which label got the full string. IF axis EQ 0 THEN BEGIN ra = value dec = 0.0 ; ticlabels,ra,1,0,ticlabs,/RA radec, ra, dec, ihr, imin, xsec, ideg, imn, xsc sd = '!Ah!N' & sm = '!Am!N' & ss = '!As!N' if (!d.name eq 'PS') and (!p.font eq 0) then begin ;Postscript fonts? sd ='!Uh!N' & sm='!Um!N' & ss='!Us!N' endif fmt = '(i2.2)' if index eq 0 then begin ticlab = string(ihr,fmt)+sd+' ' + $ string(imin,fmt) + sm +' '+ string(xsec,fmt) + ss endif else ticlab = string(xsec,'(i2.2)') + ss RETURN, ticlab ENDIF IF axis EQ 1 THEN BEGIN ra = 0.0 dec = value ; ticlabels,dec,6,0,ticlabs radec, ra, dec, ihr, imin, xsec, ideg, imn, xsc fmt = '(i2.2)' sd = "!Ao!N" & sm = "'" & ss = "''" if dec lt 0 then begin neg = -1 & sgn = '-' endif else begin neg = 1 & sgn = '' endelse if (!d.name eq 'PS') and (!p.font eq 0) then begin RtEF = '!X' sd = '!9' + string(176b) + RtEF sm = '!9' + string(162b) + RtEF ss = '!9' + string(178b) + RtEF endif if index eq 3 then begin ; this will only work if there are 7 labels ticlab = sgn+string(abs(ideg),fmt)+sd+' '+string(imn,fmt)+sm+$ ' '+string(xsc,fmt)+ss endif else ticlab = string(imn,fmt)+sm+' '+string(xsc,fmt)+ss RETURN, ticlab ; RETURN, ticlabs[0] ENDIF END