function [head, prof, dbt, jacinfo]=jacread(fin);

% function [head, prof, dbt, jacinfo]=jacread(fin);
%
% Uses rtpread to read in a finite difference jacobian RTP file,
% and re-arranges the data into a nicer format.
%
% Input:
%    fin = (string) name of input jac rtp file
%
% Output:
%    head = RTP "head" structure
%    prof = RTP "prof" structure (for one profile)
%    dbt = (njac x nlev x nchan) finite difference jacobians
%       dBT = BT_perturbed - BT_unperturbed
%    jacinfo = (njac x 2) jacobian info: (gasID, perturbation)
%

% Created 12 September 2001, Scott Hannon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

addpath /asl/matlab/h4tools

% Read input file
[head, hattr, profin, pattr] = rtpread(fin);

if (head.ptype < 1)
   disp('ERROR: input file is not a "layers" profile');
   return
end

% Determine number of profile
nlevs=profin.nlevs(1);
nlay=nlevs - 1;
nprof=length(profin.nlevs);
if (nprof < nlevs)
   disp('ERROR: input file does not seem to be a jac rtp file')
   return
end
nchan=head.nchan;


% Determine number of jacobian gases
% nprof = unperturbed + nlay*(temperature + each jac gas)
x=(nprof - 1)/nlay - 1;
ngjac=round(x);
if (abs(x - ngjac) > 0.001)
   disp('ERROR: unable to deduce number of jacobian gases')
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create output "prof" structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Location
if (isfield(profin,'plat'))
   prof.plat=profin.plat(1);
end
if (isfield(profin,'plon'))
   prof.plon=profin.plon(1);
end
if (isfield(profin,'ptime'))
   prof.ptime=profin.ptime(1);
end
if (isfield(profin,'landfrac'))
   prof.landfrac=profin.landfrac(1);
end
%
% Surface
if (isfield(profin,'stemp'))
   prof.stemp=profin.stemp(1);
end
if (isfield(profin,'spres'))
   prof.spres=profin.spres(1);
end
if (isfield(profin,'salti'))
   prof.salti=profin.salti(1);
end
%
% Emissivity & reflectivity
if (isfield(profin,'nrho'))
   prof.nrho=profin.nrho(1);
   prof.rfreq=profin.rfreq(:,1);
   prof.rho=profin.rho(:,1);
end
if (isfield(profin,'nemis'))
   prof.nemis=profin.nemis(1);
   prof.efreq=profin.efreq(:,1);
   prof.emis=profin.emis(:,1);
end
%
% Angles
if (isfield(profin,'scanang'))
   prof.scanang=profin.scanang(1);
end
if (isfield(profin,'satzen'))
   prof.satzen=profin.satzen(1);
end
if (isfield(profin,'satazi'))
   prof.satazi=profin.satazi(1);
end
if (isfield(profin,'solzen'))
   prof.solzen=profin.solzen(1);
end
if (isfield(profin,'solazi'))
   prof.solazi=profin.solazi(1);
end
%
% Profiles
if (isfield(profin,'plevs'))
   prof.plevs=profin.plevs(:,1);
end
if (isfield(profin,'plays'))
   prof.plays=profin.plays(:,1);
end
if (isfield(profin,'palts'))
   prof.palts=profin.palts(:,1);
end
prof.ptemp=profin.ptemp(:,1);
for ig=1:head.ngas
   gstr=int2str(head.glist(ig));
   eval(['prof.gas_' gstr '=profin.gas_' gstr '(:,1);']);
end
if (isfield(profin,'co2ppm'))
   prof.co2ppm=profin.co2ppm(1);
end
%
% Radiance
if (isfield(profin,'rcalc'))
   prof.rcalc=profin.rcalc(:,1);
else
   disp('ERROR: input file has no radiances')
   return
end


%%%%%%%%%%%%%%%%%%%%%%%%%%
% Deduce the jacobian info
%%%%%%%%%%%%%%%%%%%%%%%%%%
jacinfo=zeros(1 + ngjac, 2);
%
% Temperature
jacinfo(1,2)=profin.ptemp(1,2) - profin.ptemp(1,1);
%
% Amount
% Indices of all layer 1 gas amount jacobians
ind1jac=(1 + nlay + 1):nlay:nprof;
igjac=0;
% Loop over all gases
for ig=1:head.ngas
   gstr=int2str( head.glist(ig) );
   % Gas amount for layer 1
   reg1=eval(['profin.gas_' gstr '(1,1)']);
   eval(['gjunk=profin.gas_' gstr '(1,ind1jac) - reg1;']);
   ii=find(gjunk ~= 0);
   nii=length(ii);
   if (nii > 0)
      if (nii > 1)
         disp('ERROR: can not deduce amount perturbations')
         return
      else
         jacinfo(ii+1,1)=head.glist(ig);
         jacinfo(ii+1,2)=100*gjunk(ii)/reg1;
      end
   end
end


%%%%%%%%%%
% Calc dbt
%%%%%%%%%%
dbt=zeros(1+ngjac,nlay,nchan);
f=head.vchan;
bt=ones(nlay,1)*radtot(f,prof.rcalc/1000)';
%
% Temperature
ind=2:(nlay+1);
dbt(1,:,:)=radtot(f,profin.rcalc(:,ind)/1000)' - bt;
%
% Amount
for ij=1:ngjac
   ind=(ij*nlay + 2):((ij+1)*nlay + 1);
   dbt(ij+1,:,:)=radtot(f,profin.rcalc(:,ind)/1000)' -bt;
end


%%% end of file %%%
