pro contour_simple_2d_dhc_v1 ; modified from Mike Kuhn's code (contour_simple_2d) with the following changes: ; ; handling the new format of log files under xspec v12 ; allowing for a taustar vs. umax plot ; ; 1, 2, and 4 sigma confidence levels for delta c stat with ; 2 degrees of freedom sigmalevels=[2.3,4.61,6.17] datafile = dialog_pickfile(/read, filter = '*.log') minc=0.0d print, 'C stat of best fit:' read,minc ;Gets axis labels ;print, 'Axis labels: 1=sigma,wavelength 2=v_inf,u_max' print, 'Axis labels: 1=taustar,u_max 2=v_inf,u_max' axischoice=0 read,axischoice if (axischoice eq 1) then begin ; FROM MIKE'S VERSION OF THE CODE - FOR GAUSSIANS: CENTR. VS. WIDTH ; x_axislab='Wavelength ('+STRING("305B)+')' ; y_axislab='!Ms!N!7' x_axislab=textoidl('\tau_*') y_axislab=textoidl('u_{max}') endif else begin x_axislab=textoidl('u_{max}') y_axislab=textoidl('v_{\infty}') endelse hash = '' cstat=0.0d dcstat=0.0d index_x=0 x_val=0.0d index_y=0 y_val=0.0d ;transread,1,cstat,dcstat,x_val,y_val,filename=datafile transread,1,hash,cstat,dcstat,index_x,x_val,index_y,y_val,count=count,failcount=failcount,format='(A3,2F0,I0,F0,I0,F0)',filename=datafile ;finds the number of steps of the x and y variables i=0 while (y_val[i] eq y_val[0]) do i=i+1 x_steps=i y_steps=(n_elements(cstat)+1)/i print,"X steps",x_steps print,"Y steps",y_steps x_min=min(x_val) x_max=max(x_val) y_min=min(y_val) y_max=max(y_val) x=x_min+findgen(x_steps)*(x_max-x_min)/(x_steps-1) y=y_min+findgen(y_steps)*(y_max-y_min)/(y_steps-1) ;makes data into a two dimensional array cstat=reform(cstat,x_steps,y_steps,/overwrite) ; [1 2 3] ;we must be careful about the order things are written in [6 5 4] ; [7 8 9] ; ; .log files can either begin with the x-variable ascending or descending ; we must identify which case pertains, and the load the 2-D array accordingly ; sdrawkcab=x_steps-findgen(x_steps)-1 forwards=findgen(x_steps) for i=0,y_steps-1 do begin if x_val[0] lt x_val[1] then begin ; FOR ASCENDING INDEX IN FIRST SLICE (Mike's original) if ((i mod 2) eq 1) then cstat[forwards,i]=cstat[sdrawkcab,i] endif else begin ; FOR DESCENDING INDEX IN FIRST SLICE if ((i mod 2) eq 0) then cstat[forwards,i]=cstat[sdrawkcab,i] endelse endfor ;finds point with minimum c_statistic minval=cstat[0,0] minval_x=0 minval_y=0 for i=0,x_steps-1 do begin for j=0,y_steps-1 do begin if (cstat[i,j] lt minval) then begin minval=cstat[i,j] minval_x=i minval_y=j endif endfor endfor print,'Minimum c-stat of',minval,'at ',x(minval_x),y(minval_y) ;generates delta c stat array delta=cstat-minc for mkps=0,1 do begin if(mkps eq 1) then begin psopen, "tmp" ; !p.font=0 device,/times endif contour,delta,x,y,xr=[x_min,x_max],yr=[y_min,y_max],$ levels=sigmalevels,xs=1,ys=1,$ xtitle=x_axislab,ytitle=y_axislab,chars=2,thick=3 oplot,[x(minval_x),x(minval_x)],[y(minval_y),y(minval_y)],ps=2 if(mkps eq 1) then begin psclose endif endfor end