Dear Ole, Patrick and all,
now I am trying to run the stratospheric temperature retrievals in ARTS
2.2. These retrievals were running in the version 2.1 until now. I have
adapted our scripts following the demos for this version as example but
there is an error when the qpack2 is run. However, just with the error
message (see below) is difficult to find out what is going wrong. I
attach the functions that set up the Q and the Y structures and also the
inversion routine. It seems that arts is well installed because all the
demos work. I would appreciate if anyone could help to solve this problem.
Thank you very much advance!
Regards,
Fran
??? Error using ==> arts at 67
An error occured while executing ARTS. See above.
Error in ==> arts_oem_init>init_local at 411
arts( cfile );
Error in ==> arts_oem_init at 161
[Q,R,t_field,z_field,vmr_field,wind_u,wind_v,wind_w] = init_local( Q, R );
Error in ==> qpack2 at 97
[Qm,O,R,xa] = arts_oem_init( Qm, O, workfolder );
--
*****************************************************************************
Francisco Navas Guzmán, PhD
Microwave Physics Group
Division of Atmospheric radiometry
Institute of Applied Physics
University of Bern
Address: | http://www.iap.unibe.ch
Sidlerstrasse 5 | Ph: +41-31 631 50 46
CH-3012 Bern (Switzerland) | Fax: +41-31 631 37 65
*****************************************************************************
function [Y, data] = get_Y(start_time_str,end_time_str,ch_f_binned1,ch_f_binned2,bw_line1,bw_line2,exwi)
%get_Y - set up spectrum used in the retrieval
% get_Y(start_time_str,end_time_str,ch_f_binned1,ch_f_binned2) gets the TEMPERA data
% between the dates start_time_str and end_time_str (yyyy-mm-dd hh:mm:ss or
% yyyy-mm-dd) from mysql using the function get_level1_fft_from_mysq.m.
% Here the data is organized and the unwanted parts of the spectra are removed.
% In addition the measurement noise is determined from the spectra and the
% reference points for the hydrostatic equilibrium are taken from
% radiosonde measurements.
% Corinne Straub 2013-07-18
disp(sprintf('loading measurement between %s and %s',start_time_str, end_time_str))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initialize Y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y=qp2_y;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- O2 frequencies [GHz]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O2_f1=52.5424*1e9;
O2_f2=53.0669*1e9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- read level1 data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data=get_level1_fft_from_mysql(start_time_str,end_time_str);
if ~isempty(data.fft1)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- merge the spectra (including ch and f) around the two lines in one vector
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data.ch = [ch_f_binned1(1,:) ch_f_binned2(1,:)];
data.f=[ch_f_binned1(2,:) ch_f_binned2(2,:)]-45e3; % correction by Corinne Straub
data.y = [data.fft1 data.fft2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the channels around the linecenter (+-1MHz) which are
% affected by the zeeman effect
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inds_lc=find((data.f<O2_f1+1e6&data.f>O2_f1-1e6)|( data.f<O2_f2+1e6&data.f>O2_f2-1e6));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the channels outside of the frequency range used in the
% retrieval
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inds_orr=find((data.f<O2_f1-bw_line1/2)|(data.f>O2_f1+bw_line1/2)&(data.f<O2_f2-bw_line2/2)|(data.f>O2_f2+bw_line2/2));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find faulty channels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
faulty_chs=[7666:7676 8192 23254:5:23279 24579 25090:25110];
% faulty_chs=[7666:7676 8192 23254:5:23279 24579]; %to avoid gap in simulations
for k =1:length(faulty_chs)
faulty_inds(k)=find(data.ch==faulty_chs(k));
end
inds_to_remove=[inds_lc inds_orr faulty_inds];
% inds_to_remove=[inds_orr faulty_inds]; %to avoid gap in simulations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% remove 'no-good' channels
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data.y(inds_to_remove) =[];
data.f(inds_to_remove) = [];
data.ch(inds_to_remove) = [];
else
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- if there is no data in the recuired time interval return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('no level1 data found, inversion aborted!')
data=[];
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine sigma by removing a linear fit from a small area of the
% spectrum in the unbinned area but not to close to the peak
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ind= find(ch_f_binned1(1,:)>=7650 & ch_f_binned1(1,:)<=7750);
[p,s,mu]=polyfit(1:length(ind),data.fft1(ind),1);
level1_sigma=data.fft1(ind)-polyval(p,1:length(ind),[],mu);
data.sigma=std(level1_sigma);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine reference point for hydrostatic equilibrium from radiosonde
% data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
radiosonde = readpayerne(datestr(data.time,31),'morning');
if isempty(radiosonde.time)
radiosonde = readpayerne(datestr(data.time,31),'evening');
end
[val ind]=min(abs(radiosonde.pressure-100));
data.hse_p = radiosonde.pressure(ind);
data.hse_z = radiosonde.altitude(ind);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% tropospheric correction (subfunction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data.y_before_trop_corr=data.y;
data=tropospheric_correction(data,exwi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% fill the measurement vector Y (subfunction)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y = fill_Y(Y,data);
function Y = fill_Y(Y,data)
%fill_Y - fill the measurement vector used in the retrieval
% fill_Y(Y, data) - fills the structure Y with the measured data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the date
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.YEAR = year(mean(data.time));
Y.MONTH = month(mean(data.time));
Y.DAY = day(mean(data.time));
Y.HOUR = hour(mean(data.time));
Y.MINUTE = minute(mean(data.time));
Y.SECOND = second(mean(data.time));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lat/lon
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.LATITUDE = mean(data.lat);
Y.LONGITUDE = mean(data.lon);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The platform altitude is assumed to be at the ground. The zenith angle is
% set to 0 and later adjusted in the Q-structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.Z_PLATFORM = 12000; % Platform altitude
Y.ZA = 0; % Zenith angle of the measurement does not matter...
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The reference point for HSE is taken from radiosonde measurements
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.HSE_P=data.hse_p*100; % hPa->Pa
Y.HSE_Z=data.hse_z;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Spectrum with frequency axis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.F = data.f';
Y.Y = data.y';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% noise on spectrum with correction defined by Oliver Stähli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Y.TNOISE = data.sigma;
function data=tropospheric_correction(data,exwi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the indices of the first line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ind1=find(data.f<52.8e9);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine the measured tropospheric temperature for the first line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(ind1)
f_trop_corr1=[data.f_trop_corr1 data.f_trop_corr2];
Tb_trop_corr1=[data.Tb_trop_corr1 data.Tb_trop_corr2];
[p,s,mu]=polyfit(f_trop_corr1*1e9,Tb_trop_corr1,1 );
T_meas1=polyval(p,data.f(ind1),[],mu);
else
T_meas1=[];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% find the indices of the second line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ind2=find(data.f>52.8e9);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% determine the measured tropospheric temperature for the second line
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(ind2)
f_trop_corr2=[data.f_trop_corr3 data.f_trop_corr4];
Tb_trop_corr2=[data.Tb_trop_corr3 data.Tb_trop_corr4];
[p,s,mu]=polyfit(f_trop_corr2*1e9,Tb_trop_corr2,1 );
T_meas2=polyval(p,data.f(ind2),[],mu);
else
T_meas2=[];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% estimate the tropospheric temperature from the ground temperature
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
T_trop=exwi.temperature-12+273.15;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% do the tropospheric correction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t1= -([T_meas1; T_meas2]-T_trop)./(T_trop-2.7);
data.y=(data.y-T_trop*(1-t1))./t1;
function inversion(start_time_str,end_time_str)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define directories
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
top_dir = '/home/tempera/data_processing_v1/retrieval/Inversion_FFT_v4_arts22/retrieval_data/';
save_dir = '/home/tempera/data_processing_v1/retrieval/Inversion_FFT_v4_arts22/profiles/v1/';
binning_dir=('/home/tempera/data_processing_v1/calibration/fft/binning/');
addpath('/home/tempera/data_processing_v1/retrieval/matlab/');
unix(['mkdir ' save_dir]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% O2 linecenters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O2_f1=52.5424e9;
O2_f2=53.0669e9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load the binning files, convert frequency of ch_f_binned to Hz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load([binning_dir 'binning']);
load([binning_dir 'ch_f_binned']);
ch_f_binned_1(2,:)=ch_f_binned_1(2,:)*1e9;
ch_f_binned_2(2,:)=ch_f_binned_2(2,:)*1e9;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% write time in required format
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[start_time_str,start_time] = time_str_to_matlab_time(start_time_str);
[end_time_str,stop_time] = time_str_to_matlab_time(end_time_str);
end_time=start_time;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% some settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bw_line1=0e6; %bandwidth
bw_line2=240e6;
exwi=find_ground_pressure_exwi(datestr(start_time,31), datestr(end_time,31));
while end_time<stop_time
start_time_str=datestr(start_time,31);
if start_time>datenum(2012,06,04,04,00,00) && start_time<datenum(2012,06,29,12,10,00)
[start_time end_time sigma] = get_time_interval_sigma_trop_corr(start_time_str,2,exwi);
else
[start_time end_time sigma] = get_time_interval_sigma_trop_corr(start_time_str,2,exwi);
end
if datenum(start_time_str)==datenum(end_time_str)
disp('no new data')
return
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Get ground pressure from ExWi meteo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
exwi=find_ground_pressure_exwi(datestr(start_time,31), datestr(end_time,31));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Get spectrum from mysql and store it in format required by qpack
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[Y, data] = get_Y(datestr(start_time,31),datestr(end_time,31),ch_f_binned_1,ch_f_binned_2,bw_line1,bw_line2,exwi);
if ~isempty(Y.Y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Fill the Q-structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = setup_Q(top_dir,Y,data,exwi);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check that all frequencies are OK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~qp2_check_f( Q, Y, 1e3 );
error( 'Some mismatch between Q.F_BACKEND and frequencies of spectra.' );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% OEM variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O = qp2_l2( Q ); % This ensures that OEM returns the varibles needed
% to fill the L2 structure, as defined in Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define if retrieval is linear or non linear and give
% iteration method
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O.linear = false;
if ~O.linear
O.itermethod = 'GN';
O.maxiter = 10;
O.ga_start = 100;
O.ga_max =1e7;
O.ga_factor_ok = 10;
O.ga_factor_not_ok = 10;
O.stop_dx = 0.001;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Make inversion
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L2 = qpack2( Q, O, Y );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Make simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%spectra_simulation_v2(Q, O, Y);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Add some aditional parameters to the output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
L2.no_of_spectra_averaged = data.no_of_spectra_averaged;
L2.sigma=Y.TNOISE;
L2.sigma_uncorrected=data.sigma;
L2.start_time=datestr(start_time,31);
L2.end_time=datestr(end_time,31);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%- Save retrievals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~isempty(L2)
filename=sprintf('%s/retrieval_%i%.02i%.02i%.02i%.02i.mat',save_dir,L2.year,L2.month,L2.day,L2.hour,L2.minute);
save(filename,'L2');
insert_profile_into_mysql(L2,4,1);
end
end
start_time=end_time;
end
function [Q] = setup_Q(top_dir,Y,msm,exwi)
%setup_Q - set up the Q structure for the retrieval with qpack2
% setup_Q(top_dir,Y,msm,exwi) sets up the Q-structure for the retrieval
% with arts/qpack2. a documentation of the parameters needed is given in
% /opt/arts2/atmlab/retrieval/qpack2/qpack2.pdf (or similar)
% Corinne Straub 2013-07-18
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Atmlab settings
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
atmlab( 'VERBOSITY',3);
arts_xmldata_path = atmlab( 'ARTS_XMLDATA_PATH' );
arts_includes = atmlab( 'ARTS_INCLUDES' );
if isnan( arts_xmldata_path )
error( 'You need to set ARTS_XMLDATA_PATH to run this file.' );
end
if isnan( arts_includes )
error( 'You need to ARTS_INCLUDES to run this file.' );
end
%
fascod = fullfile( arts_xmldata_path, 'planets', 'Earth', 'Fascod' );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Q
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q = qarts;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% General
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.INCLUDES = { fullfile( 'ARTS_INCLUDES', 'general.arts' ), ...
fullfile( 'ARTS_INCLUDES', 'agendas.arts' ), ...
fullfile( 'ARTS_INCLUDES', 'continua.arts' ), ...
fullfile( 'ARTS_INCLUDES', 'planet_earth.arts' ) };
Q.ATMOSPHERE_DIM = 1;
Q.STOKES_DIM = 1;
Q.J_DO = true;
Q.CLOUDBOX_DO = false;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define agendas (new in Arts 2.2)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Q.WSMS_AT_START{1} = 'Copy(iy_main_agenda,iy_main_agenda__Emission)';
% Q.WSMS_AT_START{2} = 'Copy(ppath_agenda,ppath_agenda__FollowSensorLosPath)';
% Q.WSMS_AT_START{3} = 'Copy(ppath_step_agenda,ppath_step_agenda__GeometricPath)';
% Q.WSMS_AT_START{4} = 'Copy(blackbody_radiation_agenda,blackbody_radiation_agenda__Planck)';
% Q.WSMS_AT_START{5} = 'Copy(iy_space_agenda,iy_space_agenda__CosmicBackground)';
% Q.WSMS_AT_START{6} = 'Copy(iy_surface_agenda,iy_surface_agenda__UseSurfaceRtprop)';
Q.PPATH_AGENDA = { 'ppath_agenda__FollowSensorLosPath' };
Q.PPATH_STEP_AGENDA = { 'ppath_step_agenda__GeometricPath' };
Q.BLACKBODY_RADIATION_AGENDA = { 'blackbody_radiation_agenda__Planck' };
Q.IY_SPACE_AGENDA = { 'iy_space_agenda__CosmicBackground' };
Q.IY_SURFACE_AGENDA = { 'iy_surface_agenda__UseSurfaceRtprop' };
Q.IY_MAIN_AGENDA = { 'iy_main_agenda__Emission' };
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Radiative transfer
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.IY_UNIT = 'RJBT';
Q.YCALC_WSMS = { 'yCalc' };
Q.PPATH_LMAX = 2000;
%Q.PPATH_STEP_AGENDA = {};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Surface
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.Z_SURFACE = msm.alt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Absorption
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.ABS_LINES = [top_dir,'/Hitran08_o2.txt'];
Q.ABS_LINES_FORMAT = 'Hitran';
Q.ABS_LINESHAPE = 'Voigt_Drayson';
Q.ABS_LINESHAPE_FACTOR = 'quadratic';
Q.ABS_LINESHAPE_CUTOFF = -1;
Q.ABSORPTION = 'OnTheFly';
Q.ABS_NLS = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pressure grid - needs to extend below surface, at the surface we use the
% surface pressure from exwi_meteo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z_toa = 70e3;
step = 2000;
Q.P_GRID = z2p_simple( Q.Z_SURFACE-2*step : step : z_toa )';
if ~isnan(exwi.pressure)
Q.P_GRID(2)=exwi.pressure*100;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Frequency grid of forward model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.F_GRID=[Y.F(1)-35e4;Y.F;Y.F(end)+35e4];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sensor variables
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
H = qartsSensor;
H.SENSOR_NORM = true;
H.BACKEND_DO = true;
H.MIXER_DO = false;
H.F_BACKEND = Y.F;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Channel response of backend - for the FFT spectrometer we use a sinc with
% a FWHM of 32 kHz which is realistic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ch_resp = read_datafile([top_dir 'channelresponse.aa'],'MATRIX');%FFT_ACQ1_channelresponse.aa
B.name = 'Spectrometer channel response function';
B.gridnames = { 'Frequency' };
B.grids = {ch_resp(:,1)}';
B.dataname = 'Response';
B.data = ch_resp(:,2)';
H.BACKEND_CHANNEL_RESPONSE{1} = B;
clear B
%- Correlation of thermal noise (not correlated)
Q.TNOISE_C = covmat1d_from_cfun(H.F_BACKEND,[],'drc');
Q.SENSOR_DO = true;
Q.SENSOR_RESPONSE = H;
Q.ANTENNA_DIM = 1;
clear H f
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% set the azimuthal and zenith angle grid for each measurement block. The grid
% is given as an angular off-set with respect to the angles in *sensor_los*
% (Y.ZA?).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.MBLOCK_AA_GRID = {};
Q.MBLOCK_ZA_GRID = 90-msm.el;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define L2 structure (beside retrieval quantities)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.L2_EXTRA = { 'cost', 'dx', 'xa', 'y', 'yf', 'bl', 'ptz', 'mresp', ...
'e', 'eo', 'es', 'date','A','tnoise'};
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Temperature Retrieval and HSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.T.RETRIEVE = true;
Q.T.L2 = true;
Q.T.GRIDS = { Q.P_GRID(3:end), [], [] };
Q.T.ATMDATA = gf_load( fullfile(top_dir,'T_apriori') );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% apriori covariance similar to the one defined by Oliver Stähli
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cov1 =fliplr([1013e02 260e2 70e2 0.1e2 0.01e2;...
0.0025 0.04 2 2 0.0025]);
Q.T.SX = covmat1d_from_cfun( Q.T.GRIDS{1}, cov1','exp', 0.2, 0.55, @log10 );
Q.T.HSE = 'on';
Q.T.METHOD = 'analytical';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine altitudes through HSE
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.HSE.ON = true;
Q.HSE.P = 1013;
Q.HSE.ACCURACY = 0.1;
if ~Q.HSE.ON
Q.Z_ATMDATA = gf_artsxml( fullfile( arts_xmldata_path, 'climatology', ...
'cira86', 'cira86.z.xml' ), 'Altitudes', 'z_field' );
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Absorption in Troposphere and Stratosphere
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Oxygen
Q.ABS_SPECIES(1).TAG = {'O2'};
Q.ABS_SPECIES(1).RETRIEVE = false;
% Water
Q.ABS_SPECIES(2).TAG = {'H2O-PWR98'};
Q.ABS_SPECIES(2).RETRIEVE = false;
% Nitrogen
Q.ABS_SPECIES(3).TAG = {'N2-SelfContStandardType'};
Q.ABS_SPECIES(3).RETRIEVE = false;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Use Standardatmosphere separate for summer and winter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Y.MONTH>4 && Y.MONTH<10
season = 'summer';
else
season = 'winter';
end
Q.ABS_SPECIES(1).ATMDATA = gf_artsxml(fullfile( arts_xmldata_path, ...
sprintf('planets/Earth/Fascod/midlatitude-%s',season), ...
sprintf('midlatitude-%s.O2.xml',season)),'O2','vmr_field');
Q.ABS_SPECIES(2).ATMDATA = gf_load( fullfile(top_dir,'H2O_apriori') );
Q.ABS_SPECIES(3).ATMDATA = gf_artsxml(fullfile( arts_xmldata_path, ...
sprintf('planets/Earth/Fascod/midlatitude-%s',season), ...
sprintf('midlatitude-%s.N2.xml',season)),'N2','vmr_field');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Supplement tropospheric H2O with radiosonde data
% assume exponential decay from ground value for tropospheric vmr
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ind1=find(Q.ABS_SPECIES(2).ATMDATA.GRID1/100>30);
% z=p2z_simple(Q.ABS_SPECIES(2).ATMDATA.GRID1(ind1));
% vmr_trop = exwi.h2o_vmr./exp(z/2000);
%
% Q.ABS_SPECIES(2).ATMDATA.DATA(ind1)=vmr_trop;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Polyfit - a polynomial of order 3 is used for "baseline fit".
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.POLYFIT.RETRIEVE = true;
Q.POLYFIT.ORDER = 3;
Q.POLYFIT.L2 = true;
Q.POLYFIT.SX0 = 3^2;
Q.POLYFIT.SX1 = 1^2;
Q.POLYFIT.SX2 = 0.5^2;
Q.POLYFIT.SX3 = 0.1^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sinefit - there are two periods fitted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Q.SINEFIT.RETRIEVE = true;
Q.SINEFIT.PERIODS = [56.5e6; 40.9e6];
Q.SINEFIT.SX1 = 1;
Q.SINEFIT.SX2 = 1;
Q.SINEFIT.L2 = true;
_______________________________________________
qpack mailing list
[email protected]
https://www.sat.ltu.se/mailman/listinfo/qpack
