; ROC and STONE curve calculations for the IMPTAM-GOES comparison.

close,1
device,decompose=0

monthmin=fltarr(11,12)
monthmax=fltarr(11,12)
monthdays=fltarr(11,12)

first_year=2015.

; first day of values for each month: start at year 2007
monthmin(0,*)=[1,13,1,10,1,1,2,9,1,1,1,1]  ; 2007
monthmin(1,*)=[1,6,1,11,3,16,1,4,8,1,1,1]  ; 2008
monthmin(2,*)=[1,1,1,1,1,1,1,1,1,1,1,1]  ; 2009
monthmin(3,*)=[1,1,1,1,1,1,1,1,1,1,1,1]  ; 2010
monthmin(4,*)=[1,1,1,1,5,1,1,1,9,1,1,1]  ; 2011
monthmin(5,*)=[9,1,1,1,1,1,2,1,1,1,1,1]  ; 2012
monthmin(6,*)=[1,1,1,1,9,1,1,1,1,1,1,1]  ; 2013
monthmin(7,*)=[1,1,1,1,1,1,1,1,1,1,4,1]  ; 2014
monthmin(8,*)=[1,1,1,1,1,1,1,1,1,1,1,1]  ; 2015
monthmin(9,*)=[1,1,1,1,1,1,1,1,1,1,1,1]  ; 2016
monthmin(10,*)=[1,1,1,1,1,1,1,1,1,1,1,1] ; 2017

; last day of values for each month
monthmax(0,*)=[ 0,28,26,30,31,13,31,31,30,31,30,31]  ; 2007
monthmax(1,*)=[31,29,30,30,17,30,31,21,30,31,30,31]  ; 2008
monthmax(2,*)=[31,28,31,30,31,30,31,31,30,31,30,31]  ; 2009
monthmax(3,*)=[31,28,31,30,31,30,31,31,30,26,26,31]  ; 2010
monthmax(4,*)=[31,28,31,30,31,18,31,31,30,31,30,31]  ; 2011
monthmax(5,*)=[31,29,31,30,31,30,31,31,30,31,30,31]  ; 2012
monthmax(6,*)=[31,28,31,30,31,29,31,31,30,31,30,16]  ; 2013
monthmax(7,*)=[31,28,31,30,31,30,20,31,30,10,30,31]  ; 2014
monthmax(8,*)=[31,28,31,30,31,30,31,31,30,31,30,31]  ; 2015
monthmax(9,*)=[31,29,31,30,31,30,31,31,30,31,30,31]  ; 2016
monthmax(10,*)=[31,28,31,30,31,30,31,31,30,31,30,31]  ; 2017

; actual number of days in each month
monthdays(0,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2007
monthdays(1,*)= [31,29,31,30,31,30,31,31,30,31,30,31]  ; 2008
monthdays(2,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2009
monthdays(3,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2010
monthdays(4,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2011
monthdays(5,*)= [31,29,31,30,31,30,31,31,30,31,30,31]  ; 2012
monthdays(6,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2013
monthdays(7,*)= [31,28,31,30,31,30,31,31,30,31,30,31]  ; 2014
monthdays(8,*)= [ 0, 0, 0, 0, 0,30,31,31,30,31,30,31]  ; 2015
monthdays(9,*)= [31,29,31,30,31,30,31,31,30,31,30,31]  ; 2016
monthdays(10,*)=[31,28,31,30,31,30, 0, 0, 0, 0, 0, 0]  ; 2017

yeardays=total(monthdays,2)

imax_all=long(2)*long(12)*long(31)*long(24)*long(60) ; x5 = buffer! Every 5 minutes
imax_hr=long(13)*long(12)*long(31)*long(24)
imax_d=long(13)*long(12)*long(31)

;***************************************
i_read_data = 1   ; 0 don't read, 1 read

i_doSTONEstats=1    ; =1 to do the STONE curve statistics, otherwise skip just to make the plot
i_doROCstats=1    ; =1 to do the ROC curve statistics, otherwise skip just to make the plot
i_ROCnum = 2      ; number of ROC curves to calculate
     
   ; Do the stats on the 40 keV channel, in a particular MLT range

;MLT_min = 0.0
;MLT_max = 24.0

MLT_min = 3.0
MLT_max = 9.0

If MLT_min lt 10. then c_MLT_min = '0'+STRING(FORMAT='(I1)',MLT_min) $
else c_MLT_min = STRING(FORMAT='(I2)',MLT_min)
If MLT_max lt 10. then c_MLT_max = '0'+STRING(FORMAT='(I1)',MLT_max) $
else c_MLT_max = STRING(FORMAT='(I2)',MLT_max)

p_MLT_min='_MLTmin'+c_MLT_min
p_MLT_max='_MLTmax'+c_MLT_max

; Define critical threshold values

STONE_max = 6.3
STONE_min = 3.3
STONE_range = STONE_max - STONE_min
STONE_nlev = 50
STONE_levels = STONE_min + findgen(STONE_nlev+1)*STONE_range/STONE_nlev

If STONE_max lt 10. then p_STONE_max = '_STONEmax0'+STRING(FORMAT='(F4.2)',STONE_max) $
else p_STONE_max = '_STONEmax'+STRING(FORMAT='(F5.2)',STONE_max)
If abs(STONE_min) lt 10. then p_STONE_min = '_STONEmin0'+STRING(FORMAT='(F4.2)',abs(STONE_min)) $
else p_STONE_min = '_STONEmin'+STRING(FORMAT='(F5.2)',abs(STONE_min))

fileout='IMPTAM_GOES_'+p_MLT_max+P_MLT_min+'_STONE_'+p_STONE_max+p_STONE_min

ROC_max = STONE_max
ROC_min = STONE_min
ROC_range = ROC_max - ROC_min
ROC_nlev = STONE_nlev
ROC_levels = STONE_levels
data_threshold=[alog10(2.E5), alog10(3.e4)]

If ROC_max lt 10. then p_ROC_max = '_ROCmax0'+STRING(FORMAT='(F4.2)',ROC_max) $
else p_ROC_max = '_ROCmax'+STRING(FORMAT='(F5.2)',ROC_max)
If abs(ROC_min) lt 10. then p_ROC_min = '_ROCmin0'+STRING(FORMAT='(F4.2)',abs(ROC_min)) $
else p_ROC_min = '_ROCmin'+STRING(FORMAT='(F5.2)',abs(ROC_min))

fileout=fileout+'_ROC_'+p_ROC_max+p_ROC_min

              
If i_read_data eq 1 then begin

;--------------------------------
; read in IMPTAM and GOES values

time_year=intarr(imax_all)
time_month=intarr(imax_all)
time_day=intarr(imax_all) 
time_hour=intarr(imax_all) 
time_minute=intarr(imax_all)
IMPTAM_40=fltarr(imax_all)
IMPTAM_75=fltarr(imax_all)
IMPTAM_150=fltarr(imax_all)
G13_MAGED_40=fltarr(imax_all) 
G13_MAGED_75=fltarr(imax_all)  
G13_MAGED_150=fltarr(imax_all)
MLT=fltarr(imax_all)
Bz_imf=fltarr(imax_all)
By_imf=fltarr(imax_all)
Bmag_imf=fltarr(imax_all)
Vsw=fltarr(imax_all)
psw=fltarr(imax_all)
Kp=fltarr(imax_all) 
AE=fltarr(imax_all) 
SYMH=fltarr(imax_all) 

time_decimal=fltarr(imax_all)


igarb=1. ; header lines in the data file
garb=''
num_read=fltarr(20)
iread=long(-1)
iread_hr=long(-1)
iread_d=long(-1)
i_new = 0

Openr,1,'Combined_IMPTAM_MAGED_OMNIDATA.dat'
For i=1,igarb Do Readf,1,garb
while not(EOF(1)) do begin
  iread=iread+long(1)
  Readf,1,num_read
  time_year(iread)=num_read(0)
  time_month(iread)=num_read(1)
  time_day(iread)=num_read(2)
  time_hour(iread)=num_read(3)
  time_minute(iread)=num_read(4)
  IMPTAM_40(iread)=num_read(5)
  IMPTAM_75(iread)=num_read(6)
  IMPTAM_150(iread)=num_read(7)
  G13_MAGED_40(iread)=num_read(8) 
  G13_MAGED_75(iread)=num_read(9)  
  G13_MAGED_150(iread)=num_read(10)
  MLT(iread)=num_read(11)
  Bz_imf(iread)=num_read(12)
  By_imf(iread)=num_read(13)
  Bmag_imf(iread)=num_read(14)
  Vsw(iread)=num_read(15)
  psw(iread)=num_read(16)
  Kp(iread)=num_read(17) 
  AE(iread)=num_read(18) 
  SYMH(iread)=num_read(19) 

  iyear = time_year(iread) - 2007. ; so that =0 is year 2007.
  day_of_year = time_day(iread) ; January 1 is day_of_year=1
  If time_month(iread) gt 1 then day_of_year = day_of_year + total(monthdays(iyear,time_month(iread)-1))
  time_of_day = time_minute(iread) + time_hour(iread)*60 ; in minutes
  time_decimal(iread) = time_year(iread) + (day_of_year - 1)/yeardays(iyear) + time_of_day/24./60.

;  If iread mod 100 eq 0 then print,iread,time_decimal(iread),time_minute(iread),IMPTAM_40(iread),SYMH(iread)

endwhile
close,1

time_year=time_year(0:iread)
time_month=time_month(0:iread)
time_day=time_day(0:iread)
time_hour=time_hour(0:iread) 
time_minute=time_minute(0:iread)
IMPTAM_40=IMPTAM_40(0:iread)
IMPTAM_75=IMPTAM_75(0:iread)
IMPTAM_150=IMPTAM_150(0:iread)
G13_MAGED_40=G13_MAGED_40(0:iread)
G13_MAGED_75=G13_MAGED_75(0:iread)
G13_MAGED_150=G13_MAGED_150(0:iread)
MLT=MLT(0:iread)
Bz_imf=Bz_imf(0:iread)
By_imf=Bz_imf(0:iread)
Bmag_imf=Bmag_imf(0:iread)
Vsw=Vsw(0:iread)
psw=psw(0:iread)
Kp=Kp(0:iread) 
AE=AE(0:iread) 
SYMH=SYMH(0:iread)
time_decimal=time_decimal(0:iread)

IMPTAM_40_log=alog10(IMPTAM_40)
IMPTAM_75_log=alog10(IMPTAM_75)
IMPTAM_150_log=alog10(IMPTAM_150)
G13_MAGED_40_log=alog10(G13_MAGED_40)
G13_MAGED_75_log=alog10(G13_MAGED_75)
G13_MAGED_150_log=alog10(G13_MAGED_150)


endif  ; end of i_read_data check

;print,max(imptam_40),max(G13_MAGED_40),min(imptam_40),min(G13_MAGED_40),max(imptam_150),max(G13_MAGED_150),min(imptam_150),min(G13_MAGED_150)
  

;********************************************

; search through the two arrays for contingency table yes and no values

; In _ROCnSTONE version of this code, do both calculations and overplot!

If i_doSTONEstats eq 1 then begin

TruePos_STONE=fltarr(STONE_nlev+1)
FalsePos_STONE=fltarr(STONE_nlev+1)
i_good_crit_STONE=intarr(STONE_nlev+1)
Thresholds=fltarr(STONE_nlev+1)

For crit_loop=0,STONE_nlev do begin  ; begin critical value loop

critical_value = STONE_levels(crit_loop)
thresholds(crit_loop) = critical_value

crit_IMPTAM = critical_value
crit_GOES = critical_value ; <= dynamic in the STONE version

hits=long(0)
misses=long(0)
false_a=long(0)
corr_neg=long(0)

good_GOES=fltarr(iread+1)  ; here through the next for loop: _hr instead of _d
good_IMPTAM=fltarr(iread+1)
igood=long(-1)
diffsq=0.

For i_t1=0,iread do begin
  bad_value = 0
  If G13_MAGED_40(i_t1) - G13_MAGED_40(i_t1) ne 0 then bad_value = 1 ; toss out NaNs
  If IMPTAM_40(i_t1) - IMPTAM_40(i_t1) ne 0 then bad_value = 1   
  If MLT(i_t1) lt MLT_min then bad_value = 1  ; limit MLT range
  If MLT(i_t1) gt MLT_max then bad_value = 1
  If bad_value eq 0 then begin  ; exclude NaN values and focus on selected MLT range
      igood=igood+long(1)
      good_GOES(igood)=G13_MAGED_40_log(i_t1)
      good_IMPTAM(igood)=IMPTAM_40_log(i_t1)
      If G13_MAGED_40_log(i_t1) ge crit_GOES then GOES_is_crit=1 else GOES_is_crit = 0
      If IMPTAM_40_log(i_t1) ge crit_IMPTAM then IMPTAM_is_crit=1 else IMPTAM_is_crit = 0
      If GOES_is_crit eq 1 and IMPTAM_is_crit eq 1 then hits=hits + long(1)
      If GOES_is_crit eq 1 and IMPTAM_is_crit eq 0 then misses = misses + long(1)
      If GOES_is_crit eq 0 and IMPTAM_is_crit eq 1 then false_a = false_a + long(1)
      If GOES_is_crit eq 0 and IMPTAM_is_crit eq 0 then corr_neg = corr_neg + long(1)
      diffsq=diffsq+(good_GOES(igood) - good_IMPTAM(igood))^2
  endif
endfor

good_GOES=good_GOES(0:igood)
good_IMPTAM=good_IMPTAM(0:igood)



True_Pos = float(hits) / float(hits+misses)
False_Pos = float(false_a) / float(false_a+corr_neg)

TruePos_STONE(crit_loop)=True_pos
FalsePos_STONE(crit_loop)=False_pos

i_good_crit_STONE(crit_loop)=igood

;If crit_loop mod 10 eq 0 then print,gap_loop,crit_loop,critical_value,igood,True_Pos,False_Pos,hits,misses,false_a,corr_neg

endfor  ; end critical value loop

print,'   Done with that round:',i_good_crit_STONE(STONE_nlev/2),TruePos_STONE(STONE_nlev/2),FalsePos_STONE(STONE_nlev/2)

endif   ; i_doSTONEstats check


;********************************************

; search through the two arrays for contingency table yes and no values

; In _ROCnSTONE version of this code, do both calculations and overplot!

If i_doROCstats eq 1 then begin

; Do the ROC stats on the 40 keV channel, in a particular MLT range

TruePos_ROC=fltarr(ROC_nlev+1,i_ROCnum)
FalsePos_ROC=fltarr(ROC_nlev+1,i_ROCnum)
i_good_crit_ROC=intarr(ROC_nlev+1,i_ROCnum)
Thresholds=fltarr(ROC_nlev+1,i_ROCnum)


For crit_loop=0,ROC_nlev do begin  ; begin critical value loop

For ROC_loop=0,i_ROCnum-1 do begin  ; begin loop over ROC curves to calculate
  
critical_value = ROC_levels(crit_loop)
thresholds(crit_loop) = critical_value

crit_IMPTAM = critical_value
crit_GOES = data_threshold(ROC_loop)   ; static in the ROC version

hits=long(0)
misses=long(0)
false_a=long(0)
corr_neg=long(0)

good_GOES=fltarr(iread+1)  ; here through the next for loop: _hr instead of _d
good_IMPTAM=fltarr(iread+1)
igood=long(-1)
diffsq=0.

For i_t1=0,iread do begin
  bad_value = 0
  If G13_MAGED_40(i_t1) - G13_MAGED_40(i_t1) ne 0 then bad_value = 1 ; toss out NaNs
  If IMPTAM_40(i_t1) - IMPTAM_40(i_t1) ne 0 then bad_value = 1
  If MLT(i_t1) lt MLT_min then bad_value = 1  ; limit MLT range
  If MLT(i_t1) gt MLT_max then bad_value = 1
 If bad_value eq 0 then begin  ; exclude NaN values and focus on selected MLT range
      igood=igood+long(1)
      good_GOES(igood)=G13_MAGED_40_log(i_t1)
      good_IMPTAM(igood)=IMPTAM_40_log(i_t1)
      If G13_MAGED_40_log(i_t1) ge crit_GOES then GOES_is_crit=1 else GOES_is_crit = 0
      If IMPTAM_40_log(i_t1) ge crit_IMPTAM then IMPTAM_is_crit=1 else IMPTAM_is_crit = 0
      If GOES_is_crit eq 1 and IMPTAM_is_crit eq 1 then hits=hits + long(1)
      If GOES_is_crit eq 1 and IMPTAM_is_crit eq 0 then misses = misses + long(1)
      If GOES_is_crit eq 0 and IMPTAM_is_crit eq 1 then false_a = false_a + long(1)
      If GOES_is_crit eq 0 and IMPTAM_is_crit eq 0 then corr_neg = corr_neg + long(1)
      diffsq=diffsq+(good_GOES(igood) - good_IMPTAM(igood))^2
  endif
endfor

good_GOES=good_GOES(0:igood)
good_IMPTAM=good_IMPTAM(0:igood)



True_Pos = float(hits) / float(hits+misses)
False_Pos = float(false_a) / float(false_a+corr_neg)

TruePos_ROC(crit_loop,ROC_loop)=True_pos
FalsePos_ROC(crit_loop,ROC_loop)=False_pos

i_good_crit_ROC(crit_loop,ROC_loop)=igood

;If crit_loop mod 10 eq 0 then print,gap_loop,crit_loop,critical_value,igood,True_Pos,False_Pos,hits,misses,false_a,corr_neg

endfor  ; end ROC curve calc loop

endfor  ; end critical value loop

print,'   Done with that round:',i_good_crit_ROC(ROC_nlev/2,0),TruePos_ROC(ROC_nlev/2,0),FalsePos_ROC(ROC_nlev/2,0)

endif   ; i_doROCstats check

;********************************************
; Make subsampled arrays for symbol overlay of the lines

TruePos_ROC_few=fltarr(ROC_nlev+1,i_ROCnum)
FalsePos_ROC_few=fltarr(ROC_nlev+1,i_ROCnum)
TruePos_STONE_few=fltarr(ROC_nlev+1)
FalsePos_STONE_few=fltarr(ROC_nlev+1)
icount=-1
For i=0,ROC_nlev,5 do begin  ; skip however many you want, or even set it to 1
  icount=icount+1
  for ROC_loop=0,i_ROCnum-1 do begin
    TruePos_ROC_few(icount,ROC_loop)=TruePos_ROC(i,ROC_loop)
    FalsePos_ROC_few(icount,ROC_loop)=FalsePos_ROC(i,ROC_loop)
  endfor
  TruePos_STONE_few(icount)=TruePos_STONE(i)
  FalsePos_STONE_few(icount)=FalsePos_STONE(i)
endfor
TruePos_ROC_few=TruePos_ROC_few(0:icount,*)
FalsePos_ROC_few=FalsePos_ROC_few(0:icount,*)
TruePos_STONE_few=TruePos_STONE_few(0:icount)
FalsePos_STONE_few=FalsePos_STONE_few(0:icount)
  

;********************************************

printfile=1  ; set to 0 if in IDL demo mode

If printfile eq 1 then begin

  ;Print out all of these numbers to a file

  openw,1,fileout+'.dat'
  printf,1,'Filename = '+fileout+'.dat'
  printf,1,'Statistics of IMPTAM real-time results and GOES-13 values'
  printf,1,'Data threshold for ROC calculation (log flux): ',data_threshold(*)
  printf,1,''
  printf,1,format='(7A11,A20)','Threshold ',' POD-STONE ', ' POFD-STONE ', 'POD-ROC ','POFD-ROC '
  for crit_loop=0,STONE_nlev do begin
     printf,1,format='(F11.2,7F11.6)',STONE_levels(crit_loop),TruePos_STONE(crit_loop),FalsePos_STONE(crit_loop), $
        TruePos_ROC(crit_loop,*),FalsePos_ROC(crit_loop,*)
  endfor
  close,1

Endif  ; printfile check


;************************************************
; make a plot of the receiver-operator characteristic curve and scatter plot



makeplot=1
If makeplot eq 1 then begin

  !P.multi=[0,0,0]

  nplots2= 1 ; number of columns
  xgap=5.
  xleft=45.  ;xgap
  xright=15.
  xsize=275.
  xsz=xleft+xright+nplots2*(xgap+xsize)
  nplots=1   ; number of rows
  ygap=5.
  ysize=150.
  yupper=25.
  ylower=40.
  ysz=ylower+yupper+nplots*(ygap+ysize)
  y1=ysize/ysz
  xs1=xgap/xsz
  xs3=xleft/xsz
  xs2=xsize/xsz
  label_x=xs2*.1
  ys1=ygap/ysz
  ys3=ylower/ysz
  ys2=ysize/ysz
  label_y=ys2*.05
  n1=nplots*nplots2-1
  ctk1=0.8  ;.99
  ctk2=1.
  ctk3=.7
  !P.charsize=ctk1
  Legendlab=['!6STONE curve','!6ROC Curve (high thr.)','!6ROC Curve (low thr.)']


; Plot values

For ic=1,2 Do Begin
If ic Eq 1 Then Begin
 set_plot,'X'
 window,0,title=fileout 
 lth=2
Endif Else Begin
 set_plot,'ps'
 device,file=fileout+'.eps',xsize=.01*xsz,ysize=.01*ysz,/inch, $
  /port,yoff=1.,scale_factor=1.0,color=8
 lth=3
Endelse

LOADCT,39
ia=!D.table_size-1
!P.background=ia
erase,ia

;lstr=[0,1,2,3,4,5,0]
lstr=[0,0,0,0,0,0,0]
;colors=[black,red,blue,orange,green,purple,cyan]
lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5]


xx=xs3
yy=ys3

!x.range=[0.,1.]
!y.range=[0.,1.]
plot,FalsePos_STONE(*),TruePos_STONE(*),xtitle='!6Probability of False Detection', $
  ytitle='!6Prob. of Detection', $
  position=[xx,yy,xx+xs2,yy+ys2],charsize=1.,thick=lth, /nodata,  $
  linestyle=0,title='!6ROC & STONE, MLT='+c_MLT_min+' to MLT='+c_MLT_max,color=lcol(0),xstyle=1,ystyle=1
oplot,[0.,1.],[0.,1.],linestyle=1,color=lcol(5),thick=lth
oplot,FalsePos_STONE,TruePos_STONE,thick=lth,linestyle=0,color=lcol(1)
oplot,FalsePos_STONE_few,TruePos_STONE_few,thick=lth,psym=6,color=lcol(1)
for ROC_loop=0,i_ROCnum-1 do begin
  oplot,FalsePos_ROC(*,ROC_loop),TruePos_ROC(*,ROC_loop),thick=lth,linestyle=0,color=lcol(2+ROC_loop)
  oplot,FalsePos_ROC_few(*,ROC_loop),TruePos_ROC_few(*,ROC_loop),thick=lth,psym=6,symsize=.5,color=lcol(2+ROC_loop)
endfor


  a=.40
  b=.10
  c=.025
  e=.10
  d=.33
  f=.025
  For i=0,i_ROCnum Do oplot,[a,a+b],[d-i*e,d-i*e],linestyle=lstr(i+1),thick=lth,color=lcol(i+1)
  for i=0,i_ROCnum Do xyouts,a+b+c,d-i*e-f,Legendlab(i),color=0,charsize=ctk3



If ic NE 1 Then Begin
 device,/close
 set_plot,'X'
endif
endfor

;********************************************
; binning points for contour overlays

valmin=2.
valmax=7.

countbinnum = 50 ; sets resolution of the contour plots
val_range = valmax - valmin
binsize = val_range/float(countbinnum)  
count_boundaries = findgen(countbinnum+1)*binsize + valmin
countcenters = count_boundaries-0.5*binsize
countcenters(0)=valmin
countcenters(countbinnum)=valmax
countbins = intarr(countbinnum+1, countbinnum+1)
countbins(*,*) = 0

For i=0,igood do begin
  If good_goes(i) gt valmin then goes_bin = min(where(count_boundaries gt good_goes(i))) else goes_bin=-1
  If good_imptam(i) gt valmin then imptam_bin = min(where(count_boundaries gt good_imptam(i))) else imptam_bin=-1
  If goes_bin*imptam_bin ge 0 then countbins(imptam_bin,goes_bin) = countbins(imptam_bin,goes_bin)+1
;  If i mod 1000 eq 0 then print,i,good_goes(i),goes_bin,count_boundaries(goes_bin),good_imptam(i),imptam_bin,count_boundaries(imptam_bin),countbins(imptam_bin,goes_bin)
endfor

print,'binned counts: ',max(countbins)

binlines = [50,100,150,200,250,300]
binlinestyles = [0,0,0,0,0,0]  ;[1,2,1,2,1,2]
binlinecols = [1,6,4,3,2,5]


;************************************************
; make a plot of the scatterplot of SWMF and Kyoto

  For ic=1,2 Do Begin
    If ic Eq 1 Then Begin
      set_plot,'X'
      window,1,title=fileout
      lth=2
    Endif Else Begin
      set_plot,'ps'
      device,file=fileout+'_scatter.eps',color=8
      lth=3
    Endelse

    LOADCT,39
    ia=!D.table_size-1
    !P.background=ia
    erase,ia

    ;lstr=[0,1,2,3,4,5,0]
    lstr=[0,0,0,0,0,0,0]
    ;colors=[black,red,blue,orange,green,purple,cyan]
    lcol=[0,ia-1,ia/4,4*ia/5,3*ia/5,ia/10,2*ia/5]

    !x.range=[valmin,valmax]
    !y.range=[valmin,valmax]
    plot,good_imptam,good_goes,xtitle='!6IMPTAM Output',ytitle='!6GOES Observations', $
      psym=7, symsize=.2,title='!6Log10 Fluxes, MLT='+c_MLT_min+' to MLT='+c_MLT_max,color=0,xstyle=1,ystyle=1
    oplot,[valmin,valmax],[valmin,valmax],linestyle=2,color=lcol(5),thick=lth
    for ROC_loop=0,i_ROCnum-1 do oplot,[-1000.,1000.],[data_threshold(ROC_loop),data_threshold(ROC_loop)], $
      linestyle=2,color=lcol(2+ROC_loop),thick=lth
    contour,countbins,countcenters,countcenters,c_colors=lcol(binlinecols),xstyle=1,ystyle=1, $
      levels=binlines,c_linestyle=binlinestyles,/overplot,thick=lth


    If ic NE 1 Then Begin
      device,/close
      set_plot,'X'
    endif
  endfor



endif ; makeplot check


end