; .run ./prep_target_spectra_comparisons.pro
; warning, that takes a while...

!p.multi=[0,2,2,0,1]
!x.margin=[6,2]
!y.margin=[2,1]
!x.omargin=[2,0]
!y.omargin=[2,0]
!p.charsize=1.1
;whlambda=0 ; 70 mA set filename accordingly!
whlambda=2 ; 200 mA set filename accordingly

colarr=[30,80,150,250]

set_plot,'ps'
;device,file='PAPER_histo_by_structure_70_v3.eps',/encap,bits=8,/color,/times,xsi=20,ysi=12
device,file='PAPER_histo_by_structure_200_v3.eps',/encap,bits=8,/color,/times,xsi=20,ysi=12
; v2 is with the rebinning/smoothing, original w/o.  Barely any difference by eye, 
; which is all that counts here...

xtickcirc=['0','1.0','2.0','3.0','4.0']
;IVM, CIRCULAR
loadct,0
plot,indgen(5),indgen(5),xran=[0,400],yran=[0,0.02],/xst,/yst,xtickname=xtickcirc,$
         ytit='Fraction of Points, IVM',/nodata
xyouts,300,0.015,'(a)',size=2
loadct,13
for i=0, 3 do begin 
;   tmpm=(icircdl[whlambda,64:447,64:447])(where(ivmmask eq i))
   tmpm=(icircdl[whlambda,*,*])(where(ivmmask eq i))
   histo,tmpm*1000.,min=1,max=399,xx,yy,/noplot
   nn=n_elements(tmpm)
   oplot,xx,yy/float(nn),psym=10,col=colarr[i],thick=4
endfor

;MODEL, Circular
loadct,0
plot,indgen(5),indgen(5),xran=[0,400],yran=[0,0.02],/xst,/yst,xtickname=xtickcirc,$
         xtit='Fractional Polarization (%), Circular',ytit='Fraction of Points, MODEL',/nodata
xyouts,300,0.015,'(c)',size=2
loadct,13
for i=0, 3 do begin 
   tmpm=0.
   for j=2, 3 do begin
      tmpm=[tmpm,(mcircdl[whlambda,*,*,j])(where(mm.mask eq i))]
   endfor
   tmpm=tmpm[1:*]
   nn=n_elements(tmpm)
   histo,tmpm*1000.,xx,yy,min=1,max=399,/noplot
   ;oplot,xx,smoothe(yy/float(nn),2),psym=10,col=colarr[i],thick=4
   oplot,xx,yy/float(nn),psym=10,col=colarr[i],thick=4
endfor

xticklin=['0.0','0.5','1.0','1.5','2.0']
;'IVM, Linear,
loadct,0
plot,indgen(5),indgen(5),xran=[0,200],yran=[0,0.02],/xst,/yst,xtickname=xticklin,$
         /nodata
xyouts,150,0.015,'(b)',size=2
loadct,13
for i=0, 3 do begin 
   tmpm=(ilindl[whlambda,*,*])(where(ivmmask eq i))
   nn=n_elements(tmpm)
   histo,tmpm*1000.,xx,yy,min=1,max=199,/noplot
    oplot,xx,yy/float(nn),psym=10,col=colarr[i],thick=4
endfor

;'MODEL, Linear;
loadct,0
plot,indgen(5),indgen(5),xran=[0,200],yran=[0,0.02],/xst,/yst,xtickname=xticklin,$
         xtit='Fractional Polarization (%), Linear',/nodata
xyouts,150,0.015,'(d)',size=2
loadct,13
for i=0, 3 do begin 
   tmpm=0.
   for j=2, 3 do begin
      tmpm=[tmpm,(mlindl[whlambda,*,*,j])(where(mm.mask eq i))]
   endfor
   tmpm=tmpm[1:*]
   nn=n_elements(tmpm)
;   print,n_elements(tmpm)
   histo,tmpm*1000.,xx,yy,min=1,max=199,/noplot
   oplot,xx,yy/float(nn),psym=10,col=colarr[i],thick=4
endfor
device,/close

!p.multi=0
device,file='PAPER_IVM_10486_mask.eps',/color,bits=8,/times,/encap,xsi=18,ysi=16
loadct,13
contour,ivmmask,c_col=[30,80,150,250,0],/fill,lev=[0,1,2,3,5],/xst,/yst
device,/close

set_plot,'x'

reset
loadct,0


end