function [imHDR, im_data, g] = make_HDR_image(optic, t0)
  % [ix40, ix40_data] = make_HDR_image('ITMX', '2016-07-16-00-21-45');

  % look for files around the specified time
  file_tags = {'A', 'B', 'C', 'D', 'E'};
  num_file = numel(file_tags);
  
  im_data = zeros(480, 640, num_file);
  for n = 1:num_file
    im_data(:, :, n) = load_image(optic, t0, file_tags{n});
  end

  g = [30, 100, 300, 1000, 3000, 10000];
  imHDR = im_data(:, :, 1);
  for n = 1:(num_file - 1)
    %g(n) = get_gain(im_data(:, :, n), im_data(:, :, n + 1)) * g(n - 1);
    
    imHDR = merge_images(imHDR, im_data(:, :, n + 1), g(n + 1) / g(1));
  end
end

function im_data = load_image(optic, t0, file_tag)
  % make file tag
  switch optic
    case 'ITMX'
      file_mid = 'H1 ITMX (h1cam21)';
    case 'ITMY'
      file_mid = 'H1 ITMY (h1cam23)';
    case 'ETMX'
      file_mid = 'H1 ETMX (h1cam25)';
    case 'ETMY'
      file_mid = 'H1 ETMY (h1cam27)';
  end
  file_head = ['/ligo/data/camera/multiexp_' file_tag '_' file_mid '_'];
  file_tail = '.tiff';
  

  % make date stamp
  sec_per_day = 24 * 60 * 60;
  date_format = 'yyyy-mm-dd-hh-MM-SS';
  date0 = datenum(t0, date_format) - 10 / sec_per_day;

  % find file
  im_data = [];
  for n = 1:60
    file_date = datestr(date0 + n / sec_per_day, date_format);
    file_name = [file_head file_date file_tail];
    if exist(file_name, 'file') == 2
      %im_data = rot90(double(imread(file_name)), 2);
      im_data = flipud(double(imread(file_name)));
      break
    end
  end
  if isempty(im_data)
    disp(file_name)
    error('Image not found!')
  end
end

function imAB = merge_images(imA, imB, gainBA)
  % gainBA is roughly mean(imB(nn) ./ imA(nn))
  %   where nn are valid pixels
  %
  % it is assumed that gain > 1 (e.g., imB has higher gain
  %  or longer exposure time and imA)

  % start with imB with gain correction
  imAB = imB / gainBA;
  
  % use points from imA where imB is saturated
  nn = find(imB > 250);
  imAB(nn) = imA(nn);
end

function gain = get_gain(imA, imB)

  % vectorize data
  imA = imA(:);
  imB = imB(:);
  
  % find valid points
  nn = find(imA > 50 & imA < 250 & imB < 250 & imB > 50);
  imA = imA(nn);
  imB = imB(nn);

  % compute rough ratio (biased, but good starting point)
  gain = mean(imB ./ imA);
  
  % make better cut to reduce bias
  radAB = sqrt((imB / gain).^2 + imA.^2);
  nn = find(radAB > 100 & radAB < 130);
  imA = imA(nn);
  imB = imB(nn);

  % final gain
  gain = mean(imB ./ imA);  
end
