pro alpha_bb_AK19
;cal factors co:422, x:79.2, x channel at co with OD 619 

data = ''
shot_spacing = 1
plt = 0
Z = FINDGEN(506)*0.1089 - 5
cal_arr = fltarr(2000)

;for Alaska
Nfiles = 340
dir = strarr(340)
close, 5
openr, 5, 'C:\Users\jchurnside\Documents\Arctic 2019\results\clear file list.csv'
for file = 0, 339 do begin
  readf, 5, data
  dir[file] = data
endfor
out = 'C:\Users\jchurnside\Documents\Arctic 2019\results\alpha bb.csv'
close, 6
openw, 6, out
printf, 6, 'lon, lat, P(ice), alpha, sd, bb, sd, depol, file, block'

for file = 0, Nfiles-1 do begin
  rawdata = read_png(dir[file], /order)
  result = size(rawdata)
  Nshots = result[1]
  if dir[file] eq 'C:\raw data\AK 19\5-22\FL191430010.png' then Nshots = 1900
  latitude = fltarr(Nshots)
  longitude = fltarr(Nshots)
  alpha = fltarr(Nshots)
  alphac = fltarr(Nshots)
  alphax = fltarr(Nshots)
  bb = fltarr(Nshots) 
  bbc = fltarr(Nshots)
  bbx = fltarr(Nshots)
  depol = fltarr(Nshots)
  Pfit = fltarr(Nshots) 
  ice_arr = fltarr(Nshots)
  
  for shot = 0, Nshots-1, shot_spacing do begin
    SUB_READ_AUX_AK19, RAWDATA, SHOT, TIME, LAT, LON, ALT, CO_GAIN, X_GAIN, TILT, temp, inst_temp, V0, V1, ice
    ice = 0
    if min(rawdata[shot,2000:3999]) eq 0 then ice = 1
    if dir[file] eq 'C:\raw data\AK 19\5-26\FL191470032.png' and shot eq 1684 then ice = 1
    ice_arr[shot] = float(ice)
    if lon eq 0 OR ice eq 1 then continue
    latitude[shot] = lat
    longitude[shot] = -lon
    SUB_CONVERT_SHOT_AK19, RAWDATA, SHOT, CO_GAIN, X_GAIN, 'CO', ICATH
    SUB_SURFACE_CORRECT_AK19, 'sea', 0, ICATH, I2, SURFI
    if surfi eq 0 then continue             ; surface in first 500 shots. Probably clouds
    bco = 422*(I2-mean(ICATH[1900:1999])) * (300/alt)^2
    SUB_CONVERT_SHOT_AK19, RAWDATA, SHOT, CO_GAIN, X_GAIN, 'X', ICATH
    SUB_SURFACE_CORRECT_AK19, 'sea', 0, ICATH, I2, SURFI
    if surfi eq 0 then continue
    bx = 79.2*(I2-mean(ICATH[1900:1999])) * (300/alt)^2              ; without attenuator
    bt = bco+bx
    
    ; total
    logbt = alog(bt)
    ix1 = IZ(3)
    ix2 = IZ(5)      ; originally 10, then 5
    if min(bt[ix1:ix2]) lt 0 then continue
    result = linfit(Z[ix1:ix2], logbt[ix1:ix2], Yfit=Yfit, prob=prob)
    alpha[shot] = -0.5*result[1]
    bb[shot] = 2*!PI*exp(result[0])
    cal_arr[shot] = mean(bx[64:138]/bco[64:138])
    Pfit[shot] = 1-prob
    
    ; co
    logbco = alog(bco)
    result = linfit(Z[ix1:ix2], logbco[ix1:ix2], Yfit=Yfitc, prob=prob)
    alphac[shot] = -0.5*result[1]
    bbc[shot] = 2*!PI*exp(result[0])
    
    ; x
    logbx = alog(bx)
    result = linfit(Z[ix1:ix2], logbx[ix1:ix2], Yfit=Yfitx, prob=prob)
    alphax[shot] = -0.5*result[1]
    bbx[shot] = 2*!PI*exp(result[0])
    if bbx[shot] gt 0 then depol[shot] = bbx[shot]/bbc[shot]
    
   if alpha[shot] lt 0.05 then continue
   ;if alpha[shot] lt 0.05 then begin
   if alpha[shot]/bb[shot] gt 200 then begin
      PLOT, Z[0:505], bco[0:505]+1e-6, /ylog, yrange=[1e-7,1e-1], background=-1, color=0, title = shot
      OPLOT, Z[0:505], bx[0:505]+1e-6, color=255
      oplot, Z[0:505], bt[0:505]+1e-6, color=256*255
      oplot, Z[ix1:ix2], exp(Yfit)+1e-6, color=128, thick=2
      ; plot, Z[0:505], bx[0:505]/bco[0:505]
      WAIT, 0.1
      print, dir[file], shot
      stop
    endif
    
  endfor      ; end of shot loop
    
  Nblocks = fix(float(Nshots)/500)
  for block = 0, Nblocks-1 do begin
    good = where(latitude[500*block:500*block+499] gt 0, Ngood)
    if Ngood gt 5 then printf, 6, mean(longitude[500*block+good],/NAN), mean(latitude[500*block+good],/NAN), $
      total(ice_arr[500*block:500*block+499])/500.0, mean(alpha[500*block+good],/NAN), stddev(alpha[500*block+good],/NAN), $
      mean(bb[500*block+good],/NAN), stddev(bb[500*block+good],/NAN), mean(depol[500*block+good],/NAN), dir[file], block, $
      format='(8(F12.4,","),A,",",I2)'
  endfor
  
endfor
close, /all
print, mean(cal_arr), stddev(cal_arr)
end

function IZ, Z
IZ = fix((Z+5)/0.1089 + 0.5)
return, IZ
end