% TEST_ARTS_JACOBIANS
%
%    Runs a set of tests to check the implementation of Jacobians.
%
%    No options here, all changes must be done in code.
%
% FORMAT   test_arts_jacobians

% 2005-06-16   Created by Patrick Eriksson


function test_arts_jacobians


%= Select cases to run
%
% Case 1 is always included
%
do_tests = [ 1:4 ];


%= Make quite
%
verbosity = atmlab( 'VERBOSITY' );
artsverb  = atmlab( 'FMODEL_VERBOSITY' );
%
atmlab( 'VERBOSITY', 0 );
atmlab( 'FMODEL_VERBOSITY', 0 );


%= Init
%
addpath_ami;
%
Q = set_q;



%= Test 1: Just spectrum 
%
Q.JACOBIAN_DO = 0;
%
if 1
  %
  y0 = arts_y( Q );
  %
  t = [  y0(1)   min(y0)   max(y0)   y0(end)];
  e = [12.6083   11.7278   33.0063   14.7772];
  %
  if max(abs( t - e )) > 1e-3
    disp( ' ' )
    disp( 'Expected values:' )
    disp( e )
    disp( 'Obtained values:' )
    disp( t )
    error('Incorrect spectrum values obtained')
  else
    disp( 'Spectrum values with sensor : OK' );
  end
  %
end



%= Test 2: Gas jacobians, relative unit, with sensor
%
% Compare analytical and perturbation calculations
%
Q.JACOBIAN_DO               = 1;
Q.GAS_SPECIES_JAC{1}.DO     = 1;
Q.GAS_SPECIES_JAC{1}.DX     = 0.01e-6;
%
if any( do_tests == 2 )
  %
  [y,dy,Ja] = arts_y( Q );
  %
  if max(abs( y - y0 )) > 1e-4
    error( 'Spectrum vector changed when activating analytical gas jacobians');
  end
  %
  Q.GAS_SPECIES_JAC{1}.METHOD = 'perturbation';
  [y,dy,Jp] = arts_y( Q );
  %
  if max(abs( y - y0 )) > 1e-4
    error( 'Spectrum vector changed when activating analytical gas jacobians');
  end
  %
  cmpr_J( Ja, Jp, 'Gas jacobians, relative unit' );
  %
end



%= Turn off sensor
%
Q.SENSOR_RESPONSE = NaN;



%= Test 3: Gas jacobians, VMR unit
%
% Compare analytical and perturbation calculations
%
Q.GAS_SPECIES_JAC{1}.UNIT   = 'vmr';
Q.GAS_SPECIES_JAC{1}.METHOD = 'analytical';
%
if any( do_tests == 3 )
  %
  [y,dy,Ja] = arts_y( Q );
  %
  Q.GAS_SPECIES_JAC{1}.METHOD = 'perturbation';
  %
  [y,dy,Jp] = arts_y( Q );
  %
  cmpr_J( Ja, Jp, 'Gas jacobians, VMR unit' );
  %
end



%= Test 4: Gas jacobians, ND unit
%
% Compare analytical and perturbation calculations
%
Q.GAS_SPECIES_JAC{1}.UNIT   = 'nd';
Q.GAS_SPECIES_JAC{1}.METHOD = 'analytical';
%
if any( do_tests == 4 )
  %
  [y,dy,Ja] = arts_y( Q );
  %
  Q.GAS_SPECIES_JAC{1}.METHOD = 'perturbation';
  %
  [y,dy,Jp] = arts_y( Q );
  %
  cmpr_J( Ja, Jp, 'Gas jacobians, ND unit' );
  %
end




%= Reset verbosities
%
atmlab( 'VERBOSITY', verbosity );
atmlab( 'FMODEL_VERBOSITY', artsverb );


return




%----------------------------------------------------------------------------

function   cmpr_J( Ja, Jp, type );
  %
  d = full( max(max(abs( Jp-Ja ))) / max(max(abs(Ja))) ); 
  fprintf( 'Rel. difference between anal. and pert. jacobians : %.3f\n', d );
  if  d > 0.02
    keyboard
    error( sprintf('Inconsistency found for %s', type ) );
  else
    fprintf( '%s : OK\n', type );
  end
  %
return


%----------------------------------------------------------------------------

function Q = set_q


%= Init Q structures
%
Q  = qarts;
Q1 = qarts1;
H  = qarts_sensor;



%= Atmosphere
%
Q.ATMOSPHERE_DIM      = 2;
Q.USE_RAW_ATMOSPHERE  = 1;
if Q.ATMOSPHERE_DIM > 1
  Q.RAW_ATM_EXPAND_1D = 1;
end
arts_xmldata_path     = atmlab( 'ARTS_XMLDATA_PATH' );
if isnan( arts_xmldata_path )
  error('You need to ARTS_XMLDATA_PATH to run this example.');
end
Q.RAW_ATMOSPHERE      = fullfile( arts_xmldata_path, 'atmosphere', ...
                                                        'fascod', 'tropical' );
Q.P_GRID              = z2p_simple( 0e3:1e3:100e3 );
if Q.ATMOSPHERE_DIM >= 2
  Q.LAT_GRID          = -10:10;  
end


%= Spectroscopy
%
Q.GAS_SPECIES{1}{1}         = 'O3';
Q.GAS_SPECIES_JAC{1}.DO     = 0;
Q.GAS_SPECIES_JAC{1}.METHOD = 'analytical';
Q.GAS_SPECIES_JAC{1}.DX     = 0.01;
Q.GAS_SPECIES_JAC{1}.UNIT   = 'rel';
Q.GAS_SPECIES_JAC{1}.P_GRID = z2p_simple( [20e3:4e3:100e3] );
if Q.ATMOSPHERE_DIM >= 2
  Q.GAS_SPECIES_JAC{1}.LAT_GRID = -10:2:10;  
end
%
Q.GAS_SPECIES{2}{1}         = 'ClO';
Q.GAS_SPECIES{3}{1}         = 'N2O';
Q.GAS_SPECIES{4}{1}         = 'H2O';
Q.GAS_SPECIES{4}{2}         = 'H2O-MPM89';
Q.GAS_SPECIES{5}{1}         = 'N2-SelfContStandardType';
%
Q.F_GRID            = linspace( 501.18e9, 501.58e9, 201 );
Q.STOKES_DIM        = 1;
%
Q1.LINEFORMAT       = 'Arts';
Q1.LINEDATA         = fullfile( atmlab_example_data , 'lines501.4' );


%= Create absorption lookup table by ARTS1
%
Q.GAS_ABS_LOOKUP = arts_abstable_from_arts1( Q, Q1, linspace(-40,40,3) );


%= RTE 
%
Q.Z_SURFACE      = 500;
Q.SURFACE_PROP_AGENDA = ...
   { 'InterpAtmFieldToRteGps(surface_skin_t,t_field){}', ...
     'surfaceFlat{"water-liebe93"}' };
%
Q.PPATH_STEP_AGENDA = { 'ppath_stepGeometric{50e3}' };
Q.Y_UNIT            = 'RJ';
%
zplat        = 600e3;
Q.SENSOR_POS = Q.R_GEOID + zplat;
if Q.ATMOSPHERE_DIM >= 2
  Q.SENSOR_POS = [ Q.SENSOR_POS -23 ];
end
Q.SENSOR_LOS = geomztan2za( Q.R_GEOID, zplat, 25e3 );


%= Sensor
%
H.MBLOCK_ZA_GRID  = [-0.04:0.01:0.04];
H.ANTENNA_LOS     = 0;
H.ANTENNA_DIAGRAM = fullfile( atmlab_example_data, 'antenna.xml' );
Q.SENSOR_RESPONSE = H;