% Just a script (not a function) to make the final, desired
% analysis plots from analyze_narrow_features
%
% 2008/07/02 IJC @ UCLA: Making revisions for the ApJ Letter
% 2007/09/11 IJC @ UCLA

colors = 'brgcmkbrgcmkbrgcmkbrgcmkbrgcmkbrgcmk';
markers = '.ox+*sdv^ph<>.ox+*sdv^ph<>.ox+*sdv^ph<>';

if ~multiWLtest
  % Make a raw time-series plot
  tf = figure; ta = get(tf, 'children'); 
  ordToPlot = min(3, numel(output));
  lamInd = intersect([2 7], 1:numel(output{ordToPlot}.lambdas));
  figures(tf).data.x = output{ordToPlot}.hjd-2453000;
  figures(tf).data.y = [];
  for ii=1:numel(lamInd);
	normdat = output{ordToPlot}.spectraIntensity(lamInd(ii),:)./repmat(median(output{ordToPlot}.spectraIntensity(lamInd(ii),:),2),1,numel(output{ordToPlot}.hjd));
    plot(output{ordToPlot}.hjd-2453000, normdat,'.')
    hold on
    figures(tf).data.y = [figures(tf).data.y; normdat];
  end
  markPlot(tf)
  xlabel('time (HJD - 2453000)')
  ylabel('(relative flux density)')
  title('Raw extracted time series for some indices')
  yl = ylim;
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  
  
  tf = figure; ta = get(tf, 'children');
  lamInd = 7;
  normdat_det = output{ordToPlot}.detRelSpec(lamInd,:)./repmat(median(output{ordToPlot}.detRelSpec(lamInd,:),2),1,numel(output{ordToPlot}.hjd));
  trend_in_data = normdat - normdat_det;
  
  plot(output{ordToPlot}.hjd-2453000, trend_in_data,'o', output{ordToPlot}.hjd-2453000, trend, '--')
  xlabel('time (HJD - 2453000)')
  ylabel('(relative flux density)')
  title('Timeseries, wavelength-normalized trend')
  ylim(yl);
  figures(tf).data.x = output{ordToPlot}.hjd-2453000;
  figures(tf).data.y = normdat_det;
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  
  tf = figure; ta = get(tf, 'children');
  plot(output{ordToPlot}.hjd-2453000, normdat_det,'.')
  xlabel('time (HJD - 2453000)')
  ylabel('(relative flux density)')
  title('Detrended timeseries, wavelength-normalized, trend removed')
  ylim(yl);
  figures(tf).data.x = output{ordToPlot}.hjd-2453000;
  figures(tf).data.y = normdat_det;
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
end


for jj=1:numel(output)
  tf = figure; ta = get(tf, 'children');
  hold on
  leg = '';
  for ii=1:numel(output{jj})
    errorbar(lambdas{jj}, output{jj}(ii).normDiff, output{jj}(ii).normErr,'.')
    leg = strvcat(leg, num2str(output{jj}(ii).refWL));
  end
  markPlot(tf);
  errorbar(lambdas{jj}, spectralTemplate{jj}, templateErrors{jj}, 'b-', 'linewidth', 2)
  leg = strvcat(leg, 'Spectral Template');
  plot(ll{jj}, modelTemplate{jj}, 'r-o', 'linewidth', 2)
  leg = strvcat(leg, 'Model Template');
  tl = legend(leg);
  ylabel('( \Delta N )')
  xlabel('Wavelength (\mum)')
  grid on
  figures(tf).data.x{1} = lambdas{jj};
  figures(tf).data.x{2} = ll{jj};
  figures(tf).data.y{1} = spectralTemplate{jj};
  figures(tf).data.y{2} = modelTemplate{jj};
  figures(tf).data.err{1} = templateErrors{jj};
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  figures(tf).legend = get(tl, 'string');
 
end

if ~multiWLtest
  tf = figure; ta = get(tf, 'children'); 
  hold on
  for jj=1:numel(output)
    leg = '';
    errorbar(lambdas{jj}, output{jj}.normDiff, output{jj}.normErr, 'bo', 'linewidth', 2)
    errorbar(lambdas{jj}, output{jj}.normDiffSys, output{jj}.normErrSys, 'go', 'linewidth', 2)
    leg = strvcat(leg, ['Ord. ' num2str(orders(jj)) ': ref. \lambda' ...
			' = ' num2str(output{jj}.refWL)], ...
		  ['Ord. ' num2str(orders(jj)) ': systematics']);
    figures(tf).data.x{jj} = lambdas{jj}(:)';
    figures(tf).data.y{jj} = [output{jj}.normDiff(:)'; output{jj}.normDiffSys(:)'];
    figures(tf).data.err{jj} = [output{jj}.normErr(:)'; output{jj}.normErrSys(:)'];
  end
  leg = strvcat(leg, 'Spectral Template');
  ylabel('( \Delta N )')
  xlabel('Wavelength (\mum)')
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  figures(tf).legend = get(tl, 'string');
end


if multiWLtest
  ordToPlot = intersect(orders, 4:5);
  for kk=1:numel(ordToPlot)
    tf = figure; ta = get(tf, 'children');
    jj = find(orders==ordToPlot(kk));
    switch ordToPlot(kk)
     case 4
      refWLind = 3;
     case 5
      refWLind = 4;
    end
    for ii=1:numel(output{jj})
      if ii==refWLind
	plot(lambdas{jj}, output{jj}(ii).normDiff, '-o', 'linewidth',2, ...
	     'color', [1 .9025 0])
	hold on
      else
	plot(lambdas{jj}, output{jj}(ii).normDiff, ['-o' colors(ii)])
	hold on
      end
      figures(tf).data.x{ii} = lambdas{jj};
      figures(tf).data.y{ii} = output{jj}(ii).normDiff;
    end
    plot(lambdas{jj}, output{jj}(refWLind).normDiff, '-o', 'linewidth',2, ...
	 'color', [1 .9025 0])
    tl = legend(num2str(lambdas{jj}'));
    xlabel('Wavelength (\mum)')
    ylabel('(\Delta N)')
    figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
    figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
    figures(tf).labels.title = get(get(ta, 'title'), 'string');
    figures(tf).legend = get(tl, 'string');
  end
else

% $$$   if stitchTogether
% $$$     tf = figure; ta = get(tf, 'children');
% $$$     hold on
% $$$     xlim([1.9 4.2])
% $$$     ylim([-2 3]*1e-3);                        
% $$$     leg = '';
% $$$     for ii=1:numel(orders)             
% $$$       errorbar(lambdas{ii}, spectralTemplate{ii}, templateErrors{ii}, 'b-')
% $$$       leg = strvcat(leg, ['SpeX Order #' num2str(ii)]);
% $$$     end
% $$$     for ii=1:numel(orders)             
% $$$       plot(ll{ii}, modelTemplate{ii}, 'r--')
% $$$       leg = strvcat(leg, ['Model, Order #' num2str(ii)]);
% $$$     end
% $$$     markPlot(tf)
% $$$     
% $$$     xlabel('Wavelength [um]')
% $$$     ylabel('Normalized A-B (arbitrary units)')
% $$$     grid on
% $$$     title('"Stitched" synthetic relative normalized difference spectrum')
% $$$     tl = legend(leg);
% $$$   end
  
  [ll2, out] = model_norm_relspec(inOptions,'--wlrange=[2.1 8]');
  tf = figure; ta = get(tf, 'children');
  plot(ll2, out.starSpectra./(pi*starRad^2), 'b-',ll2, out.planetSpectra/(pi*planetRad^2),'-.r')
  figures(tf).data.x{1} = ll2;
  figures(tf).data.y{1} = [out.starSpectra(:)'./(pi*starRad^2); out.planetSpectra(:)'/(pi*planetRad^2)];
  
  if synthSpec
    hold on
    for ii=1:numel(output)
      errorbar(lambdas{ii}, synthSpectrumCont{ii}/(pi*planetRad^2), synthErr{ii}/(pi*planetRad^2),'-g')
      figures(tf).data.x{ii+1} = lambdas{ii};
      figures(tf).data.y{ii+1} = synthSpectrumCont{ii}/(pi*planetRad^2);
      figures(tf).data.errs{ii+1} = synthErr{ii}/(pi*planetRad^2);
    end
  end
  xlim([1.9 4.2])
  xlabel('Wavelength [um]')
  ylabel('Flux [erg/s/Hz/cm^2]')
  tl = legend('Star Model', 'Planet Model', 'Synthesized Planet Spectrum');
  title('Model spectra used for HD189733')
  grid on
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  figures(tf).legend = get(tl, 'string');
  
  tf = figure; ta = get(tf, 'children');
  plot(ll2, 1e7*out.planetSpectra/(pi*planetRad^2),'-.r')
  figures(tf).data.x{1} = ll2;
  figures(tf).data.y{1} = out.planetSpectra(:)'/(pi*planetRad^2);
  if synthSpec
    hold on
    for ii=1:numel(output)
      errorbar(lambdas{ii}, 1e7*synthSpectrumCont{ii}/(pi*planetRad^2), 1e7*synthErr{ii}/(pi*planetRad^2),'xg')
      figures(tf).data.x{ii+1} = lambdas{ii};
      figures(tf).data.y{ii+1} = synthSpectrumCont{ii}/(pi*planetRad^2);
      figures(tf).data.errs{ii+1} = synthErr{ii}/(pi*planetRad^2);
    end
  end
  xlim([1.9 4.2])
  xlabel('Wavelength [um]')
  ylabel('Intensity [10^{-7} erg/s/Hz/cm^2]')
  grid off
  figures(tf).labels.x = get(get(ta, 'xlabel'), 'string');
  figures(tf).labels.y = get(get(ta, 'ylabel'), 'string');
  figures(tf).labels.title = get(get(ta, 'title'), 'string');
  figures(tf).legend = get(tl, 'string');
  
end