Dear all, Dear Patrick
Thank you! I didn't expect this to happen so quickly! I had the time to
play around a bit and it seems to work nicely. Attached is an extended
example cfile with 3D atmosphere, wind, and freq shift.
Two _minor_ things I noticed:
1) If the frequency grid f_grid is very narrow at certain points (~< 100
kHz, intentionally or by mistake with VectorInsertGridPoints) the
sensor.cc throws an error:
arts: /opt/arts-dev/arts/src/sensor.cc:1447: void
integration_func_by_vecmult(VectorView, ConstVectorView,
ConstVectorView, ConstVectorView): Assertion `is_increasing( x_f_in )'
failed.
Aborted (core dumped)
To reproduce, increase the number of points in line 56 from 101 to 102.
If that really is a limitation, catching that earlier (in
sensor_checkedCalc ?) would be helpful.
2) If the f_backend has more than approx. 6000 channels (but I guess
that number depends on the amount of RAM), the setup of the covariance
matrices and the retrieval eat up all memory and take very long. With
atmlab, I could use the 10k channels of our instrument directly. That
was very convenient, but I can work around that. Simon know about this
already, just writing it down again.
Best regards,
Jonas
On 24.09.2018 11:41, Patrick Eriksson wrote:
Hi all, Hi Jonas, Hi Richard,
All: If you do retrievals and using most recent version, you need to
do some small changes. The benefit is that you now also can retrieve
frequency corrections. See the ChangeLog message for a summary.
Jonas: With the commit I just made it should be possible to retrieve a
frequency shift using ARTS's OEM. The demo cfile should be pretty
close to what you want to do. That is, for details see:
controlfiles/artscomponents/oem/TestOEM.arts
I have done some testing, and all looked OK. Tell me if it works for
you or not.
Richard: I also made a small preparation if we ever will test to
retrieve spectroscopic variables. The mapping from the x-vector to
spectroscopic data is planned to be handled by x2artsSpectrosocopy (so
far a dummy WSM)
Bye,
Patrick
On 2018-09-18 15:18, Jonas Hagen wrote:
Dear Patrick
Thanks, of course this would be the proper approach and I hope that
you decide to implement instrument variables soon! This is the last
piece in the way of operational WIRA-C Wind retrievals with ARTS (and
thus without Matlab). Until then I will tinker with my patch.
Best regards,
Jonas
On 18.09.2018 10:27, Patrick Eriksson wrote:
Dear Jonas,
Some quick feedback, I am on a conference.
The variables that you can retrieve using ARTS-OEM so far are mainly
atmospheric quantities. To handle instrument variables I maonly left
for the future.
A constant frequency switch could be handled by shifting the
transitions as you tried to do. But frequency shifts are in fact an
instrument parameter, and moving the transitions will not work for
frequency stretch. So it seems to be time to implement a general way
to handle instrument variables. I will discuss with Simon, that is
also here at the conference.
For the moment I suggest that you do repeated linear inversions. And
adjusting your instrument frequencies after each linear inversion.
This should work with your extension of xaStandard. After some
iterations turn off the frequency switch and make a final inversion.
Kind regards,
Patrick
On 2018-09-17 17:22, Jonas Hagen wrote:
Hello ARTS Developers,
I'm trying to retrieve the Frequnecy Shift along with Wind with the
ARTS internal retrieval. To my understanding, this should work, but
support in xaStandard and x2artsStandard WSMs is missing and
results in an error: "Found a retrieval quantity that is not yet
handled by internal retrievals: Frequency"
For xaStandard(), the a priori Frequency Shift could easily be set
to zero after line 793 of m_oem.cc along with baseline and pointing.
For x2artsStandard(), maybe a new WSV would make sense (f_shift)
and the inversion_iterate_agenda would then call
abs_linesShiftFrequency(f_shift), similar to the baseline stuff?
I tried to implement it myself but got stuck with
jacobian_quantities and indices in x2artsStandard().
Best regards,
Jonas Hagen
_______________________________________________
arts_dev.mi mailing list
arts_dev...@lists.uni-hamburg.de
https://mailman.rrz.uni-hamburg.de/mailman/listinfo/arts_dev.mi
--
Jonas Hagen
Institute of Applied Physics
University of Bern
Sidlerstr. 5, CH-3012 Bern, Switzerland
Tel: +41 31 631 5045
#DEFINITIONS: -*-sh-*-
#
# Author: Patrick Eriksson
Arts2 {
INCLUDE "general/general.arts"
INCLUDE "general/continua.arts"
INCLUDE "general/agendas.arts"
INCLUDE "general/planet_earth.arts"
# Agendas to use
#
Copy( abs_xsec_agenda, abs_xsec_agenda__noCIA )
Copy( propmat_clearsky_agenda, propmat_clearsky_agenda__OnTheFly )
Copy( iy_main_agenda, iy_main_agenda__Emission )
Copy( iy_space_agenda, iy_space_agenda__CosmicBackground )
Copy( iy_surface_agenda, iy_surface_agenda__UseSurfaceRtprop )
Copy( ppath_agenda, ppath_agenda__FollowSensorLosPath )
Copy( ppath_step_agenda, ppath_step_agenda__GeometricPath )
# Basic settings
#
AtmosphereSet3D
IndexSet( stokes_dim, 1 )
# Observations
#
# Set to 1 (W), 2 (W/E) or 4 (W/E/N/S)
IndexCreate( n_obs )
IndexSet( n_obs, 2 )
# Variables to define frequency grid and
#
NumericCreate( f0 )
NumericCreate( df_start )
NumericCreate( df_end )
VectorCreate( fgrid1 )
VectorCreate( fgrid2 )
VectorCreate( fgrid3 )
IndexCreate( nf )
#
NumericSet( f0, 110.836e9 )
#NumericSet( f0, 142.17504e9 )
#
# Create broad, coarse frequency grid
NumericSet( df_start, -100e6 )
NumericSet( df_end, 100e6 )
VectorNLinSpace( fgrid1, 51, df_start, df_end )
#
# Create a narrow, fine frequency grid
NumericSet( df_start, -20e6 )
NumericSet( df_end, 20e6 )
VectorNLinSpace( fgrid2, 101, df_start, df_end )
#
# Create final f_grid
VectorInsertGridPoints( f_grid, fgrid1, fgrid2 )
VectorAddScalar( f_grid, f_grid, f0 )
#ReadXML(f_grid, "f_grid.xml")
#VectorAddScalar( f_grid, f_grid, f0 )
# Define spectrometer
#
NumericCreate( f_resolution )
NumericCreate( f_start )
NumericCreate( f_end )
#
# Use limited bandwidth, because we need high resolution for wind,
# but above approx 6k channels, memory fills up quickly.
NumericSet( f_start, -30e6 )
NumericSet( f_end, 30e6 )
NumericSet( f_resolution, 12e3 )
#
NumericAdd( f_start, f_start, f0 )
NumericAdd( f_end, f_end, f0 )
# Pressure grid
#
IndexCreate( n_p_ret )
IndexCreate( n_p )
VectorCreate( p_ret_grid )
#
IndexSet( n_p, 361 )
IndexSet( n_p_ret, 81 )
#
VectorNLogSpace( p_grid, n_p, 500e2, 0.1 )
VectorNLogSpace( p_ret_grid, n_p_ret, 500e2, 0.1 )
# Spectroscopy
#
abs_speciesSet( species=[ "O3" ] )
#
ArrayOfLineshapeSpecCreate( abs_lineshapeDefine )
abs_lineshapeDefine( abs_lineshapeDefine, "Voigt_Kuntz6", "VVH", 750e9 )
#
ReadXML( abs_lines, "testdata/ozone_line.xml" ) # 110 GHz line
#ReadXML( abs_lines, "./Perrin_O3_142.xml" ) # 142 GHz line
abs_lines_per_speciesCreateFromLines
#
IndexSet(abs_f_interp_order, 1) # no effect for OnTheFly propmat
# Latitude / Longitude
IndexCreate(n_lat)
IndexCreate(n_lon)
IndexSet(n_lon, 2)
IndexSet(n_lat, 2)
#
VectorSet( lat_true, [0] )
VectorSet( lon_true, [0] )
#
VectorNLinSpace(lat_grid, n_lat, -3, 3)
VectorNLinSpace(lon_grid, n_lon, -3, 3)
#
MatrixSetConstant( z_surface, n_lat, n_lon, 10e3 )
# Atmosphere (a priori)
#
AtmRawRead( basename = "testdata/tropical" )
AtmFieldsCalcExpand1D
#
NumericCreate(u_wind)
NumericCreate(v_wind)
NumericSet(u_wind, 60)
NumericSet(v_wind, -40)
Tensor3SetConstant(wind_u_field, n_p, n_lat, n_lon, u_wind)
Tensor3SetConstant(wind_v_field, n_p, n_lat, n_lon, v_wind)
# Apply HSE
#
NumericSet( p_hse, 100e2 )
NumericSet( z_hse_accuracy, 0.5 )
#
atmfields_checkedCalc
#
z_fieldFromHSE
# Sensor pos/los/time
#
VectorCreate(pos_alt)
VectorCreate(pos_lat)
VectorCreate(pos_lon)
VectorCreate(los_za)
VectorCreate(los_aa)
IndexCreate(tmp1)
ArrayOfIndexCreate(tmp2)
IndexStepDown(tmp1, n_obs)
ArrayOfIndexLinSpace(tmp2, 0, tmp1, 1)
VectorSetConstant(pos_alt, n_obs, 15e3)
VectorSetConstant(pos_lat, n_obs, 0)
VectorSetConstant(pos_lon, n_obs, 0)
VectorSetConstant(los_za, n_obs, 68)
VectorSet(los_aa, [90, -90, 0, 180])
Select(los_aa, los_aa, tmp2)
Delete(tmp1)
Delete(tmp2)
Matrix3ColFromVectors(sensor_pos, pos_alt, pos_lat, pos_lon)
Matrix2ColFromVectors(sensor_los, los_za, los_aa)
#
VectorSetConstant( sensor_time, n_obs, 0 )
# True f_backend
#
VectorLinSpace( f_backend, f_start, f_end, f_resolution )
nelemGet( nf, f_backend )
Print(nf)
# Assume a Gaussian channel response
#
VectorCreate( fwhm )
VectorSetConstant( fwhm, 1, f_resolution )
backend_channel_responseGaussian( fwhm = fwhm )
# With a frequency shift retrieval, we must use sensor_response_agenda
#
FlagOn( sensor_norm )
#
AgendaSet( sensor_response_agenda ){
AntennaOff
sensor_responseInit
# Among responses we only include a backend
sensor_responseBackend
}
#
AgendaExecute( sensor_response_agenda )
# RT
#
NumericSet( ppath_lmax, -1 )
StringSet( iy_unit, "RJBT" )
# Deactive parts not used (jacobian activated later)
#
jacobianOff
cloudboxOff
# Perform tests
#
abs_xsec_agenda_checkedCalc
propmat_clearsky_agenda_checkedCalc
atmfields_checkedCalc
atmgeom_checkedCalc
cloudbox_checkedCalc
sensor_checkedCalc
# Simulate "measurement vector"
#
yCalc
#
# Start on retrieval specific part
#
# Some vaiables
#
VectorCreate(vars)
SparseCreate(sparse_block)
MatrixCreate(dense_block)
# Start definition of retrieval quantities
#
retrievalDefInit
#
nelemGet( nelem, p_ret_grid )
nelemGet( nf, y )
# Add ozone as retrieval quantity
#
retrievalAddAbsSpecies(
species = "O3",
unit = "vmr",
g1 = p_ret_grid,
g2 = lat_true,
g3 = lon_true
)
VectorSetConstant( vars, nelem, 1e-12 )
DiagonalMatrix( sparse_block, vars )
covmat_sxAddBlock( block = sparse_block )
# Add zonal wind
#
retrievalAddWind(
component = "u",
g1 = p_ret_grid,
g2 = lat_true,
g3 = lon_true
)
VectorSetConstant( vars, nelem, 1e5 )
DiagonalMatrix( sparse_block, vars )
covmat_sxAddBlock( block = sparse_block )
# Add meridional wind
#
#retrievalAddWind(
# component = "v",
# g1 = p_ret_grid,
# g2 = lat_true,
# g3 = lon_true
#)
#VectorSetConstant( vars, nelem, 100 )
#DiagonalMatrix( sparse_block, vars )
#covmat_sxAddBlock( block = sparse_block )
# Add a frquency shift retrieval
#
retrievalAddFreqShift(
df = 100e3
)
VectorSetConstant( vars, 1, 1e10 )
DiagonalMatrix( sparse_block, vars )
covmat_sxAddBlock( block = sparse_block )
# Add a baseline fit
#
retrievalAddPolyfit(
poly_order = 1
)
VectorSetConstant( vars, n_obs, 5 )
DiagonalMatrix( sparse_block, vars )
covmat_sxAddBlock( block = sparse_block )
VectorSetConstant( vars, n_obs, 1 )
DiagonalMatrix( sparse_block, vars )
covmat_sxAddBlock( block = sparse_block )
# Define Se and its invers
#
VectorSetConstant( vars, nf, 1e-2 )
DiagonalMatrix( sparse_block, vars )
covmat_seAddBlock( block = sparse_block )
#
VectorSetConstant( vars, nf, 1e+2 )
DiagonalMatrix( dense_block, vars )
covmat_seAddInverseBlock( block = dense_block )
# End definition of retrieval quantities
#
retrievalDefClose
# x, jacobian and yf must be initialised (or pre-calculated as shown below)
#
VectorSet( x, [] )
VectorSet( yf, [] )
MatrixSet( jacobian, [] )
# Or to pre-set x, jacobian and yf
#
#Copy( x, xa )
#MatrixSet( jacobian, [] )
#AgendaExecute( inversion_iterate_agenda )
# Iteration agenda
#
AgendaSet( inversion_iterate_agenda ){
Ignore(inversion_iteration_counter)
# Map x to ARTS' variables
x2artsAtmAndSurf
x2artsSensor # No need to call this WSM if no sensor variables retrieved
# To be safe, rerun some checks
atmfields_checkedCalc
atmgeom_checkedCalc
# Calculate yf and Jacobian matching x.
yCalc( y=yf )
# Add baseline term (no need to call this WSM if no sensor variables
retrieved)
VectorAddVector( yf, yf, y_baseline )
# This method takes cares of some "fixes" that are needed to get the Jacobian
# right for iterative solutions. No need to call this WSM for linear
inversions.
jacobianAdjustAndTransform
}
# Let a priori be off with 0.5 ppm
#
Tensor4AddScalar( vmr_field, vmr_field, 0.5e-6 )
# Reset Wind
#
Tensor3SetConstant(wind_u_field, n_p, n_lat, n_lon, 0)
Tensor3SetConstant(wind_v_field, n_p, n_lat, n_lon, 0)
# Add a baseline
# Offset 1 K
#
VectorAddScalar( y, y, 1 )
# Introduce a frequency error
#
VectorAddScalar( f_backend, f_backend, -150e3 )
# Calculate sensor_reponse (this time with assumed f_backend)
#
AgendaExecute( sensor_response_agenda )
# Create xa
#
xaStandard
# Run OEM
OEM( method = "gn",
max_iter = 5,
display_progress = 1,
stop_dx = 0.1,
lm_ga_settings = [10,2,2,100,1,99])
#
#OEM( method = "lm",
# max_iter = 20,
# display_progress = 1,
# lm_ga_settings = [100.0, 2.0, 2.0, 10.0, 1.0, 1.0])
#
x2artsAtmAndSurf
x2artsSensor
#Print( oem_errors, 0 )
Print( x, 0 )
# Compute averaging kernel matrix
#
avkCalc
# Compute smoothing error covariance matrix
#
covmat_ssCalc
# Compute observation system error covariance matrix
#
covmat_soCalc
# Extract observation errors
#
retrievalErrorsExtract
#WriteXML( "ascii", f_backend, "f.xml" )
#WriteXML( "ascii", y, "y.xml" )
}
_______________________________________________
arts_users.mi mailing list
arts_users.mi@lists.uni-hamburg.de
https://mailman.rrz.uni-hamburg.de/mailman/listinfo/arts_users.mi
