;
pro fig4_v01
;
defsysv, '!tkref', 273.16d
defsysv, '!tkmelt', 273.15d
;
common plotting, xsize_inch, ysize_inch, width_inch, width_normal, height_inch, $
 height_normal, height_spacer_inch, height_spacer_normal, width_spacer_inch, width_spacer_normal, $
 yoffset_inch, yoffset_normal, xoffset_inch, xoffset_normal
common normal_coordinates, x_coord_min_normal,y_coord_min_normal,x_coord_max_normal,y_coord_max_normal
;
tvlct, [   0, 255,   0,   0, 255,   0, 150, 255], $
       [   0,   0,   0, 255,   0, 255, 150, 140], $
       [   0,   0, 255,   0, 255, 255, 150,  35]
black           = 0
gray            = 6
;
;
;..temp. - aerosol (Method #1) fit coefficients from Tab. 1 (revision)
;
fit_1_a    = exp(-14.89d)
fit_1_b    = 4.79d
fit_1_c    = 0.0076d
fit_1_d    = 0.86d
;
;..temp. - aerosol (D10) fit coefficients from Tab. 1
;
fit_D10_a    = exp(-9.73d)
fit_D10_b    = 3.33d
fit_D10_c    = 0.0264d
fit_D10_d    = 0.0033d
;
circle_size       = 0.02
;
lu1              = 1
infilename       = 'revised.csv'
;
;..find out how big the file is
;
openr, lu1, infilename
input_string     = ''
ndata            = 0
while not eof(lu1) do begin
 readf, lu1, input_string
 ndata = ndata + 1
endwhile
close, lu1
print, ndata
;
;..there are two header lines
;
openr, lu1, infilename
readf, lu1, input_string
print, input_string
readf, lu1, input_string
print, input_string
input_data = fltarr(9,ndata-2)
readf, lu1, input_data
close, lu1
yyyymmdd      = reform(input_data(0,*))
hhmmss0       = reform(input_data(1,*))
hhmmss1       = reform(input_data(2,*))
hhmmss2       = reform(input_data(3,*))
hhmmss3       = reform(input_data(4,*))
n05           = reform(input_data(5,*))
N             = reform(input_data(6,*))
Tlow          = reform(input_data(7,*))
tMP           = reform(input_data(8,*))
;
;..plotting
;
xrange    = [-34., -14.]
xticks    = 2
xminor    = 10
yrange    = [0.0,   18.]
yticks    = 3
yminor    = 6
;
xsize_inch            =  8.5
ysize_inch            =  11.
width_inch            =  3.0
width_normal          = width_inch / xsize_inch
height_inch           =  3.0
height_normal         = height_inch / ysize_inch
height_spacer_inch    =  1.0
height_spacer_normal  = height_spacer_inch / ysize_inch
width_spacer_inch     =  1.0
width_spacer_normal   = width_spacer_inch / xsize_inch
yoffset_inch          =  1.0
yoffset_normal        = yoffset_inch / ysize_inch
xoffset_inch          =  2.0
xoffset_normal        = xoffset_inch / xsize_inch
;
set_plot, 'ps'
!p.font    =    0
device, /times, font_size = 12, filename = 'fig4_v01.ps', $
 /inches, yoffset = 0., xoffset = 0., xsize = xsize_inch, ysize = ysize_inch, bits_per_pixel = 8, /color
;
x_coord_min_normal    = xoffset_normal
x_coord_max_normal    = x_coord_min_normal + width_normal
y_coord_min_normal    = (ysize_inch - 1.*height_inch) / ysize_inch - 0.*height_spacer_normal - yoffset_normal
y_coord_max_normal    = y_coord_min_normal + height_normal
!p.position           = [x_coord_min_normal,y_coord_min_normal,x_coord_max_normal,y_coord_max_normal]
;
xyouts, x_coord_max_normal - 0.02, y_coord_max_normal - 0.02, 'a', /normal
;
plot, /nodata, /noerase, fltarr(1), ytitle = 'D > 0.5 !9m!7m Concentration, n!D0.5!N, sccm!U-1!N', yrange=yrange, $
 xrange=xrange, xstyle=1, ystyle = 1, yticks = yticks, xticklen = -0.02, yticklen = -0.02, $
 yminor = yminor, xticks = xticks, xminor = xminor, xtitle = 'Minimum Streamline Temperature, T!Dlow!N, !Uo!NC'
box_coordinates_x = [-33., -33., -29., -29., -33.]
box_coordinates_y = [ 1.5,  3.0,  3.0,  1.5,  1.5]
polyfill, box_coordinates_x, box_coordinates_y, color=gray
make_circles, xrange, yrange, Tlow, n05, black, circle_size
;
;..put black Xs in the circles with N-measured greater than 40 sL^-1
;
index = where(N ge 40, index_count)
oplot, Tlow[index], n05[index], symsize=0.4, psym=7, thick = 1, color = black
;
;..evaluate the ratio of fitted method 1 divided by fitted p14
;
N_fit_P14  = fit_1_a * (!tkref-(Tlow[index]+!tkmelt))^fit_1_b * n05[index]^(fit_1_c*(!tkref-(Tlow[index]+!tkmelt))+fit_1_d)
N_fit_D10  = fit_D10_a * (!tkref-(Tlow[index]+!tkmelt))^fit_D10_b * n05[index]^(fit_D10_c*(!tkref-(Tlow[index]+!tkmelt))+fit_D10_d)
ratio      = N_fit_P14 / N_fit_D10
stats      = moment(ratio, sdev = my_sdev)
print, stats[0], my_sdev
print, yyyymmdd[index], format='(i15)'
;
;..checking the count
;
index = where(Tlow le xrange[1] and Tlow ge xrange[0], index_count)
print, index_count
index = where(n05 le yrange[1] and n05 ge yrange[0], index_count)
print, index_count
index = where(Tlow le xrange[1] and Tlow ge xrange[0] and $
               n05 le yrange[1] and  n05 ge yrange[0], index_count)
print, index_count
;
;..here is the selection of data from the set of 80
;..note the "lt" when comparing to Peng's analysis
;
index = where(Tlow lt max(box_coordinates_x) and Tlow ge min(box_coordinates_x) and $
               n05 lt max(box_coordinates_y) and  n05 ge min(box_coordinates_y), index_count)
print, index_count
for i = 0, index_count - 1 do begin
 print, Tlow[index[i]], n05[index[i]], N[index[i]], tMP[index[i]]
endfor
;
print, correlate(alog(tMP[index]), alog(N[index]))
result = linfit(alog(tMP[index]), alog(N[index]))
afit = result[0]
bfit = result[1]
;
x_coord_min_normal    = xoffset_normal
x_coord_max_normal    = x_coord_min_normal + width_normal
y_coord_min_normal    = (ysize_inch - 2.*height_inch) / ysize_inch - 1.*height_spacer_normal - yoffset_normal
y_coord_max_normal    = y_coord_min_normal + height_normal
!p.position           = [x_coord_min_normal,y_coord_min_normal,x_coord_max_normal,y_coord_max_normal]
;
xyouts, x_coord_max_normal - 0.02, y_coord_max_normal - 0.02, 'b', /normal
;
xrange    = [0., 400.]
xticks    = 2
xminor    = 2
yrange    = [0.0,  30.]
yticks    = 3
yminor    = 1
;
plot, /nodata, /noerase, fltarr(1), ytitle = 'N!DIC!N, sL!U-1!N', yrange=yrange, $
 xrange=xrange, xstyle=1, ystyle = 1, yticks = yticks, xticklen = -0.02, yticklen = -0.02, $
 yminor = yminor, xticks = xticks, xminor = xminor, xtitle = 'Mixed-phase Time, t!DMP!N, s'
circle_size       = 0.02
make_circles, xrange, yrange, tMP[index], N[index], black, circle_size
;
xplot = dindgen(400) + 1
yplot = exp(afit + bfit*alog(xplot))
oplot, xplot, yplot, thick = 4, color = black
;
xyouts, x_coord_min_normal + 0.02, y_coord_max_normal - 0.02, $
 'Fitted Eqn.: ln(N!DIC!N) = c!D1!N + c!D2!N' +'!9'+string(215B)+'!7'+'ln(t!DMP!N)', $
 /normal, charsize = 1.0, color = black
xyouts, x_coord_min_normal + 0.02, y_coord_max_normal - 0.04, $
 'c!D1!N = '+string(afit,format='(f5.2)'), /normal, charsize = 1.0, color = black
xyouts, x_coord_min_normal + 0.02, y_coord_max_normal - 0.06, $
 'c!D2!N = '+string(bfit,format='(f5.2)'), /normal, charsize = 1.0, color = black
xyouts, x_coord_min_normal + 0.02, y_coord_max_normal - 0.08, $
 'r = ' + string(correlate(alog(tMP[index]), alog(N[index])), format='(f5.2)') + '  p = XX', $
 /normal, charsize = 1.0, color = black
;
resetps2x
;
end
;
;
pro make_circles, xrange, yrange, x, y, color, circle_size
dx_dy    = (xrange(1) - xrange(0)) / (yrange(1) - yrange(0))
for i = 0, n_elements(x) - 1 do begin
 xcenter  = x(i)
 ycenter  = y(i)
 yheight  = (yrange(1) - yrange(0)) * circle_size
 radian   = 2.*!pi*findgen(101) / 100.
 xcircle  = (yheight/2.)*cos(radian) * dx_dy + xcenter
 ycircle  = (yheight/2.)*sin(radian) + ycenter
 oplot, xcircle, ycircle, thick = 2, color = color
 index_bias = n_elements(xcircle)/2.
; for j = 0, n_elements(xcircle)/2 - 1 do begin
;  polyfill, [xcircle(j), xcircle(j+1), xcircle(j+1), xcircle(j)], $
;             [ycircle(j),ycircle(j), ycircle(j+index_bias), ycircle(j+index_bias)], color = color
; endfor
endfor
end