Open Build Service

Changes of Revision 13

atmos.spec Changed

atmos-2.97.2.tar.gz/configure -> atmos-2.97.3.tar.gz/configure Changed

@@ -1,6 +1,6 @@
 #! /bin/sh
 # Guess values for system-dependent variables and create Makefiles.
-# Generated by GNU Autoconf 2.63 for atmos 2.97.2.
+# Generated by GNU Autoconf 2.63 for atmos 2.97.3.
 #
 # Report bugs to <mass-general@curl.sai.msu.ru>.
 #
@@ -596,8 +596,8 @@
 # Identity of this package.
 PACKAGE_NAME='atmos'
 PACKAGE_TARNAME='atmos'
-PACKAGE_VERSION='2.97.2'
-PACKAGE_STRING='atmos 2.97.2'
+PACKAGE_VERSION='2.97.3'
+PACKAGE_STRING='atmos 2.97.3'
 PACKAGE_BUGREPORT='mass-general@curl.sai.msu.ru'
 
 ac_unique_file="src/atmos.h"
@@ -1293,7 +1293,7 @@
   # Omit some internal or obsolete options to make the list less imposing.
   # This message is too long to be a string in the A/UX 3.1 sh.
   cat <<_ACEOF
-\`configure' configures atmos 2.97.2 to adapt to many kinds of systems.
+\`configure' configures atmos 2.97.3 to adapt to many kinds of systems.
 
 Usage: $0 OPTION... VAR=VALUE...
 
@@ -1359,7 +1359,7 @@
 
 if test -n "$ac_init_help"; then
   case $ac_init_help in
-     short | recursive ) echo "Configuration of atmos 2.97.2:";;
+     short | recursive ) echo "Configuration of atmos 2.97.3:";;
    esac
   cat <<\_ACEOF
 
@@ -1449,7 +1449,7 @@
 test -n "$ac_init_help" && exit $ac_status
 if $ac_init_version; then
   cat <<\_ACEOF
-atmos configure 2.97.2
+atmos configure 2.97.3
 generated by GNU Autoconf 2.63
 
 Copyright (C) 1992, 1993, 1994, 1995, 1996, 1998, 1999, 2000, 2001,
@@ -1463,7 +1463,7 @@
 This file contains any messages produced by compilers while
 running configure, to aid debugging if configure makes a mistake.
 
-It was created by atmos $as_me 2.97.2, which was
+It was created by atmos $as_me 2.97.3, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
 
   $ $0 $@
@@ -2312,7 +2312,7 @@
 
 # Define the identity of the package.
  PACKAGE='atmos'
- VERSION='2.97.2'
+ VERSION='2.97.3'
 
 
 cat >>confdefs.h <<_ACEOF
@@ -5455,7 +5455,7 @@
 # report actual input values of CONFIG_FILES etc. instead of their
 # values after options handling.
 ac_log="
-This file was extended by atmos $as_me 2.97.2, which was
+This file was extended by atmos $as_me 2.97.3, which was
 generated by GNU Autoconf 2.63.  Invocation command line was
 
   CONFIG_FILES    = $CONFIG_FILES
@@ -5518,7 +5518,7 @@
 _ACEOF
 cat >>$CONFIG_STATUS <<_ACEOF || ac_write_fail=1
 ac_cs_version="\\
-atmos config.status 2.97.2
+atmos config.status 2.97.3
 configured by $0, generated by GNU Autoconf 2.63,
   with options \\"`$as_echo "$ac_configure_args" | sed 's/^ //; s/\\""\`\$/\\\\&/g'`\\"

atmos-2.97.2.tar.gz/configure.in -> atmos-2.97.3.tar.gz/configure.in Changed

atmos-2.97.2.tar.gz/src/covariance_decomposition.h -> atmos-2.97.3.tar.gz/src/covariance_decomposition.h Changed

atmos-2.97.2.tar.gz/src/file_output.cpp -> atmos-2.97.3.tar.gz/src/file_output.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.cpp 200 2011-02-03 14:45:13Z matwey $
+    $Id: file_output.cpp 207 2011-03-01 19:02:17Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -45,17 +45,17 @@
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "M " << time << " Background measurements" << std::endl;
 }
-void file::atmosphere( const time_type& time, double free_seeing, double free_seeing_err, double dimm_seeing, double dimm_seeing_err, double dimm_turbulence, double dimm_turbulence_err, double isoplanatic_angle, double isoplanatic_angle_err, double altitude_eff, double altitude_eff_err, double desi_time_constant, double desi_time_constant_err, double time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err ) {
+void file::atmosphere( const time_type& time, double free_seeing, double free_seeing_rerr, double dimm_seeing, double dimm_seeing_rerr, double dimm_turbulence, double dimm_turbulence_rerr, double isoplanatic_angle, double isoplanatic_angle_rerr, double altitude_eff, double altitude_eff_rerr, double desi_time_constant, double desi_time_constant_rerr, double time_constant, double time_constant_rerr, double moment_0, double moment_0_rerr, double moment_2, double moment_2_rerr ) {
 	static datum skipped;
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "A " << time << " "
-		<< std::fixed << datum( free_seeing, free_seeing_err/free_seeing ) << " " << datum( dimm_seeing, dimm_seeing_err/dimm_seeing ) << " "
-		<< std::scientific << datum( moment_0, moment_0_err/moment_0 ) << " " << datum( dimm_turbulence, dimm_turbulence_err/dimm_turbulence ) << " "
-		<< std::fixed << datum( altitude_eff, altitude_eff_err/altitude_eff ) << " " << skipped << " "
-		<< datum( isoplanatic_angle, isoplanatic_angle_err/isoplanatic_angle ) << " "
-		<< std::scientific << datum( moment_2, moment_2_err/moment_2 ) << " "
-		<< std::fixed << datum( desi_time_constant, desi_time_constant_err/time_constant ) << " "
-		<< std::fixed << datum( time_constant, time_constant_err/time_constant ) << std::endl;
+		<< std::fixed << datum( free_seeing, free_seeing_rerr ) << " " << datum( dimm_seeing, dimm_seeing_rerr ) << " "
+		<< std::scientific << datum( moment_0, moment_0_rerr ) << " " << datum( dimm_turbulence, dimm_turbulence_rerr ) << " "
+		<< std::fixed << datum( altitude_eff, altitude_eff_rerr ) << " " << skipped << " "
+		<< datum( isoplanatic_angle, isoplanatic_angle_rerr ) << " "
+		<< std::scientific << datum( moment_2, moment_2_rerr ) << " "
+		<< std::fixed << datum( desi_time_constant, desi_time_constant_rerr ) << " "
+		<< std::fixed << datum( time_constant, time_constant_rerr ) << std::endl;
 }
 void file::profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const vector_type& cn2, const vector_type& ecn2, const vector_type& grid ){
 	assert( cn2.size() == ecn2.size() );

atmos-2.97.2.tar.gz/src/file_output.h -> atmos-2.97.3.tar.gz/src/file_output.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.h 200 2011-02-03 14:45:13Z matwey $
+    $Id: file_output.h 207 2011-03-01 19:02:17Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -50,7 +50,7 @@
 	void comment( const std::string& comment );
 	void normal( const time_type& time, double mz );
 	void background( const time_type& time );
-	void atmosphere( const time_type& time, double free_seeing, double free_seeing_err, double dimm_seeing, double dimm_seeing_err, double dimm_turbulence, double dimm_turbulence_err, double isoplanatic_angle, double isoplanatic_angle_err, double altitude_eff, double altitude_eff_err, double desi_time_constant, double desi_time_constant_err, double time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err );
+	void atmosphere( const time_type& time, double free_seeing, double free_seeing_rerr, double dimm_seeing, double dimm_seeing_rerr, double dimm_turbulence, double dimm_turbulence_rerr, double isoplanatic_angle, double isoplanatic_angle_rerr, double altitude_eff, double altitude_eff_rerr, double desi_time_constant, double desi_time_constant_rerr, double time_constant, double time_constant_rerr, double moment_0, double moment_0_rerr, double moment_2, double moment_2_rerr );
 	void profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const vector_type& cn2, const vector_type& ecn2, const vector_type& grid );
 	void residual( const time_type& time, const std::string& method, const vector_type& o_csi );
 	void fluxes( const time_type& time, const vector_type& result_flux );

atmos-2.97.2.tar.gz/src/indices.h -> atmos-2.97.3.tar.gz/src/indices.h Changed

atmos-2.97.2.tar.gz/src/integrals.cpp -> atmos-2.97.3.tar.gz/src/integrals.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.cpp 200 2011-02-03 14:45:13Z matwey $
+    $Id: integrals.cpp 210 2011-03-02 18:59:06Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -36,8 +36,73 @@
 namespace atmos {
 namespace integrals {
 
+value_type seeing( const value_type& moment0 ) {
+	const value_type power = 3.0/5.0;
+	const value_type k_lambda = 1.1486983549970351 * 15.848931924611138; // lambda^-1/5
+	const value_type k = 5.307 * 206264.8;
+
+	return k*k_lambda*std::pow( moment0, power );
+}
+value_type effective_altitude( const value_type& moment0, const value_type& moment1 ) {
+	const value_type k = 0.001; // km in meter
+
+	return k*moment1/moment0;
+}
+value_type isoplanatic_angle( const value_type& moment53 ) {
+	const value_type power = -3.0/5.0;
+	const value_type k_lambda = 0.43527528164806206 * 6.3095734448019363e-08; // lambda^6/5
+	const value_type k = 0.058 * 206264.8;
+
+	return k*k_lambda*std::pow( moment53, power );
+}
+value_type time_constant_desi( const value_type& desi ) {
+	const value_type power = -3.0/5.0;
+	const value_type k = 0.175;
+
+	return k*std::pow( desi, power );
+}
+value_type time_constant( const value_type& wind_moment0, const value_type& moment0 ) {
+	const value_type power = -0.1;
+	const value_type k_lambda = 0.43527528164806206 * 6.3095734448019363e-08; // lambda^6/5
+	const value_type k = 0.058 * 1e3;
+	const value_type k_wind = 6e6 / 6;
+
+	return k*k_lambda*std::pow( moment0, power )/std::sqrt( k_wind*wind_moment0 );
+}
+
+std::pair< value_type, value_type > seeing( const value_type& moment0, const value_type& moment0_rerr ) {
+	return std::make_pair( seeing(moment0), 3.0/5.0*moment0_rerr );
+}
+std::pair< value_type, value_type > effective_altitude( const value_type& moment0, const value_type& moment0_rerr, const value_type& moment1, const value_type& moment1_rerr ) {
+	return std::make_pair( effective_altitude(moment0,moment1), std::sqrt(moment0_rerr*moment0_rerr+moment1_rerr*moment1_rerr) );
+}
+std::pair< value_type, value_type > isoplanatic_angle( const value_type& moment53, const value_type& moment53_rerr ) {
+	return std::make_pair( isoplanatic_angle(moment53), 3.0/5.0*moment53_rerr );
+}
+std::pair< value_type, value_type > time_constant_desi( const value_type& desi, const value_type& desi_rerr ) {
+	return std::make_pair( time_constant_desi(desi), 3.0/5.0*desi_rerr );
+}
+std::pair< value_type, value_type > time_constant( const value_type& wind_moment0, const value_type& wind_moment0_rerr, const value_type& moment0, const value_type& moment0_rerr ) {
+	return std::make_pair( time_constant( wind_moment0, moment0), std::sqrt(0.1*moment0_rerr*moment0_rerr + 0.5*wind_moment0_rerr*wind_moment0_rerr) );
+}
+
 namespace {
 
+vector_type init_moment_powers() {
+	vector_type powers(4);
+	powers(moment_0)=0.0;
+	powers(moment_53)=5.0/3.0;
+	powers(moment_1)=1.0;
+	powers(moment_2)=2.0;
+	return powers;
+}
+
+vector_type init_wind_moment_powers() {
+	vector_type powers(1);
+	powers(moment_0)=0.0;
+	return powers;
+}
+
 // FIXME:
 vector_type reorder_flat( const atmos::indices::indices_type& indices ) {
 	assert( indices.size1() == atmos::channels );
@@ -51,70 +116,38 @@
 }
 
 }
-const double KM2M = 1000;
-const double POWSEE = 0 ;       // Power of Cn2 moment for seeing
-const double POWISP = 5.0/3.0;  // Power of Cn2 moment for isoplanatic angle
-const double POWEFF = 1.0;      // Power of Cn2 moment for Effective altitude
-const double POWISK = 2.0;      // Cn2 second moment power
-
-const double MINRELW = 6e-4 ;   // Minimal ratio of the singular value to the maximal one
-
-static const double powers = { POWSEE, POWISP, POWEFF, POWISK };
-const unsigned int npower = sizeof(powers)/sizeof(powers0);
-
-static const double wind_powers = { 0.0 };
-const unsigned int wind_npower = sizeof(wind_powers)/sizeof(wind_powers0);
-
-static const double HBOUND = 0.38;// Height of the boundary layer kilometers to fit the Cn2-integral
-static const double HTOP = 25.0; // Height of the hight layer kilometers
-
-static const int TOPLAYERS =  4; // Number of rejected high-altitude layers
-const double JPOW_SEE = 0.6 ;
-const double LPOW_SEE = -0.2;
-const double LEFF_SEE = 1.1486983549970351;
-const double KS = 1.73588e+07;  // Calibration constant in seeing formula (AstL 45, p.399):
-                                //\f$ seeing'' = KS * \lambda_{eff}^{-1/5} M_0^{3/5} *
-                                //\sec^{-1}(z)\f$, where \em lambda_eff is in mkm*/
-
-const double JPOW_ISP = -0.6;   // Power of Cn-5/3th moment for isoplan. angle
-const double LPOW_ISP = 1.2;    // Power of lambda for isoplan. angle
-const double LEFF_ISP = 0.43527528164806206;
-const double KP = 0.000756348;  // Calibration constant in isoplanatic angle formula (AstL 45, p.399)
-                                // \f$ \theta_0'' = KP * (M_{5/3} /\lambda_{eff}^2)^{-3/5} * \sec^{8/5}(z) \f$, where
-                                // \f$ KP = (2.905*(2\pi)^2*(MKM2M)^{-2})^{-3/5}*206265 \f$*/
-const double KT = 0.175;        // Atmospheric time constant empirical calibration constant for DESI for 1,3ms
-                                // integrations in diameter=2cm apertures with -1km shift
-const double KL = 0.058 * 2.7464013582652968e-08;
-
-const double SPOW_TC  = -0.6;   // Power of DESI for atm. time constant
 
-functor_factory::dcf_matrix_type functor_factory::calculate( const weight_type& wf, unsigned int npower, const double* powers ) {
+functor_factory::dcf_matrix_type functor_factory::calculate( const weight_type& wf, const vector_type& powers ) {
 	using namespace boost::numeric::ublas;
 	using namespace utils::algo;
 	using namespace utils::ublas;
 
+	const double max_cond = 1667;
+	const double bound_altitude = 0.38; // Height of the boundary layer kilometers to fit the Cn2-integral
+	const double top_altitude   = 25.0; // Height of the hight layer kilometers
+
 	const vector<double> grid( utils::ublas::vector_generator< lg_generator > ( lg_generator(0.0, 0.04, 0.25), 50 ) );
 
 	typedef dcf_matrix_type            result_type;
 	typedef dcf_matrix_type::size_type size_type;
 
-	vector<double>::const_iterator lower = std::find_if( grid.begin(), grid.end(), std::bind2nd( std::greater< double >(), HBOUND) );
-	vector<double>::const_reverse_iterator higher = std::find_if( grid.rbegin(), grid.rend(), std::bind2nd( std::less< double >(), HTOP) );
+	vector<double>::const_iterator lower = std::find_if( grid.begin(), grid.end(), std::bind2nd( std::greater< double >(), bound_altitude) );
+	vector<double>::const_reverse_iterator higher = std::find_if( grid.rbegin(), grid.rend(), std::bind2nd( std::less< double >(), top_altitude) );
 	/* FIXME: Weights::grid_type::const_iterator */
 
 	// both layer altitudes and grid are compile-time contants, so it is assertion.
 	assert( lower != grid.end() );
 	assert( higher != grid.rend() );
 
-	result_type res( npower, atmos::indices_size );
+	result_type res( powers.size(), atmos::indices_size );
 
 	vector_range< const vector<double> > h( grid, range( lower.index(), higher.index()+1 ) );
 	matrix< double > wfv = vector_of_vector_to_matrix( project( grid, wf) );
 
-	for( unsigned int ipow = 0; ipow < npower; ipow++ ) {
-		double moment_power = powersipow;
+	for( vector_type::size_type ipow = 0; ipow < powers.size(); ipow++ ) {
+		double moment_power = powers(ipow);
 		double weight_power = moment_power / 2.0;
-		double scale = pow( KM2M, moment_power );
+		double scale = pow( 1000.0, moment_power );
 
 		vector<double> alt    = utils::ublas::exp( log( h )*moment_power ) * scale;
 		vector<double> weight = utils::ublas::exp( log( h )*weight_power );
@@ -132,7 +165,7 @@
 		sd.apply( left, right );
 
 		matrix_vector_slice< matrix<double> > singular( diag( A ) ); // it is guaranteed that singular values are sorted
-		matrix_vector_slice< matrix<double> >::const_iterator small_singular = std::find_if( singular.begin(), singular.end(), std::bind2nd( std::less< double >(), singular(0)*MINRELW) );
+		matrix_vector_slice< matrix<double> >::const_iterator small_singular = std::find_if( singular.begin(), singular.end(), std::bind2nd( std::less< double >(), singular(0)/max_cond) );
 		matrix_vector_slice< matrix<double> >::const_iterator it = singular.begin();
 		alt = prod( left, alt );
 		for( ; it != small_singular; ++it ) {
@@ -162,16 +195,19 @@
 functor_factory::~functor_factory() {
 }
 const functor_factory::functor_type& functor_factory::operator() ( const std::string& spectral_type, const weight_type& null_wf, const weight_type& second_wf ) {
+	static const vector_type moment_powers = init_moment_powers();
+	static const vector_type wind_moment_powers = init_wind_moment_powers();
+
 	storage_type::const_iterator it = m_storage.find( spectral_type );
 	if( it != m_storage.end() )
 		return it->second;
 
-	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectral_type, functor_type( calculate( null_wf, npower, powers ), calculate( second_wf, wind_npower, wind_powers ) ) ) );
+	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectral_type, functor_type( calculate( null_wf, moment_powers ), moment_powers, calculate( second_wf, wind_moment_powers ), wind_moment_powers ) ) );
 
 	return ret.first->second;
 }
 
-functor::result_type functor::operator() ( const indices_type& indices, const indesis_type& indesis, const indices_type& wind_indices, double mz ) const {
+functor::result_type functor::operator() ( const indices_type& indices, const indices_type& wind_indices, const indesis_type& desi, double mz ) const {
 	typedef vector< double >       vector_type;
 	typedef vector_type::size_type size_type;
 	
@@ -183,53 +219,32 @@
 	aints = prod( m_indices_dcf, aints );
 
 	assert( npower == aints.size() );
-	for( size_type i = 0; i < npower; ++i ){
-		aints(i) = aints(i) / std::pow( mz, 1.0 + powersi );
+	for( size_type i = 0; i < m_powers.size(); ++i ){
+		aints(i) = aints(i) / std::pow( mz, 1.0 + m_powers(i) );
+		if( !std::isfinite(aints(i)) || aints(i) < 0.0 )
+			aints(i) = 0.0;
 	}
 
 	vector_type wind_aints = reorder_flat(wind_indices);
 	wind_aints = prod( m_wind_indices_dcf, wind_aints );
 		
 	assert( wind_npower == wind_aints.size() );
-	for( size_type i = 0; i < wind_npower; ++i ){
-		wind_aints(i) = wind_aints(i) / std::pow( mz, 1.0 + wind_powersi );
+	for( size_type i = 0; i < m_wind_powers.size(); ++i ){
+		wind_aints(i) = wind_aints(i) / std::pow( mz, 1.0 + m_wind_powers(i) );
+		if( !std::isfinite(wind_aints(i)) || wind_aints(i) < 0.0 )
+			wind_aints(i) = 0.0;
 	}
 
-	return result_type( KS * LEFF_SEE * std::pow( aints(0), JPOW_SEE ),
-			    KP * LEFF_ISP * std::pow( aints(1), JPOW_ISP ),
-			    0.001 * std::pow( aints(1) / aints(2), 1.0/( POWISP-POWEFF ) ),
-			    KT * std::pow( indesis(0), SPOW_TC ),
-			    1e3/*ms*/ * KL * std::pow( aints(0), -0.1 ) / std::sqrt( 6e6 / 6 * wind_aints(0) ),
-			    aints(0),
-			    aints(3) );
+	assert( desi.size() > 0 );
+
+	return boost::make_tuple( aints, wind_aints, ( !std::isfinite(desi(0)) || desi(0) < 0 ) ? 0 : desi(0) );
 }
 
 storage::storage( size_type reserve ):
 	table_storage< detail::sample >( reserve ),
-	m_free_seeing_column( *this, boost::mem_fn(&detail::sample::free_seeing) ),
-	m_isoplanatic_angle_column( *this, boost::mem_fn(&detail::sample::isoplanatic_angle) ),
-	m_altitude_eff_column( *this, boost::mem_fn(&detail::sample::altitude_eff) ),
-	m_desi_time_constant_column( *this, boost::mem_fn(&detail::sample::desi_time_constant) ),
-	m_time_constant_column( *this, boost::mem_fn(&detail::sample::time_constant) ),
-	m_moment_0_column( *this, boost::mem_fn(&detail::sample::moment_0) ),
-	m_moment_2_column( *this, boost::mem_fn(&detail::sample::moment_2) ) {
-}
-
-namespace filter {
-	
-bool regular::operator() ( const sample_type& x ) const {
-	return ! ( is_regular( x.free_seeing ) &&
-		is_regular( x.isoplanatic_angle ) &&
-		is_regular( x.altitude_eff ) &&
-		is_regular( x.desi_time_constant ) &&
-		is_regular( x.time_constant ) &&
-		is_regular( x.moment_0 ) &&
-		is_regular( x.moment_2 ) );
-}
-bool regular::is_regular( double x ) const {
-	return ! ( std::isnan( x ) || std::isinf( x ) || ( x < 0 ) );
-}
-
+	m_moments_column( *this, boost::mem_fn(&detail::sample::moments) ),
+	m_wind_moments_column( *this, boost::mem_fn(&detail::sample::wind_moments ) ),
+	m_desi_column( *this, boost::mem_fn(&detail::sample::desi) ) {
 }
 
 }}

atmos-2.97.2.tar.gz/src/integrals.h -> atmos-2.97.3.tar.gz/src/integrals.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.h 200 2011-02-03 14:45:13Z matwey $
+    $Id: integrals.h 210 2011-03-02 18:59:06Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -36,26 +36,41 @@
 	typedef vector< value_type > vector_type;
 	typedef matrix< value_type > matrix_type;
 
+	typedef vector_type moments_type;
+	typedef atmos::indices::indesis_type::value_type desi_type;
+
+	enum moment_power {
+		moment_0 = 0,
+		moment_53 = 1,
+		moment_1 = 2,
+		moment_2 = 3
+	};
+
+	value_type seeing( const value_type& moment0 );
+	value_type effective_altitude( const value_type& moment0, const value_type& moment1 );
+	value_type isoplanatic_angle( const value_type& moment53 );
+	value_type time_constant_desi( const value_type& desi );
+	value_type time_constant( const value_type& wind_moment0, const value_type& moment0 );
+
+	std::pair< value_type, value_type > seeing( const value_type& moment0, const value_type& moment0_rerr );
+	std::pair< value_type, value_type > effective_altitude( const value_type& moment0, const value_type& moment0_rerr, const value_type& moment1, const value_type& moment1_rerr );
+	std::pair< value_type, value_type > isoplanatic_angle( const value_type& moment53, const value_type& moment53_rerr );
+	std::pair< value_type, value_type > time_constant_desi( const value_type& desi, const value_type& desi_rerr );
+	std::pair< value_type, value_type > time_constant( const value_type& wind_moment0, const value_type& wind_moment0_rerr, const value_type& moment0, const value_type& moment0_rerr );
+
 namespace detail {
 	struct sample {
-		typedef atmos::integrals::value_type value_type;
+		typedef atmos::integrals::moments_type moments_type;
+		typedef atmos::integrals::desi_type desi_type;
 
-		value_type free_seeing;
-		value_type isoplanatic_angle;
-		value_type altitude_eff;
-		value_type desi_time_constant;
-		value_type time_constant;
-		value_type moment_0;
-		value_type moment_2;
-
-		sample( const value_type& free_seeing,
-			const value_type& isoplanatic_angle,
-			const value_type& altitude_eff,
-			const value_type& desi_time_constant,
-			const value_type& time_constant,
-			const value_type& moment_0,
-			const value_type& moment_2 ):
-			free_seeing(free_seeing), isoplanatic_angle(isoplanatic_angle), altitude_eff(altitude_eff), desi_time_constant(desi_time_constant), time_constant(time_constant), moment_0(moment_0), moment_2(moment_2) {
+		moments_type moments;
+		moments_type wind_moments;
+		desi_type desi;
+
+		sample( const moments_type& moments,
+			const moments_type& wind_moments,
+			const desi_type& desi ):
+			moments(moments), wind_moments(wind_moments), desi(desi) {
 		}
 	};
 }
@@ -64,51 +79,26 @@
 	public table_storage< detail::sample > {
 public:
 	typedef storage self_type;
-private:
-	typedef table_storage_column_getter< self_type, detail::sample::value_type > common_column_type;
 public:
-	typedef common_column_type free_seeing_column_type;
-	typedef common_column_type isoplanatic_angle_column_type;
-	typedef common_column_type altitude_eff_column_type;
-	typedef common_column_type desi_time_constant_column_type;
-	typedef common_column_type time_constant_column_type;
-	typedef common_column_type moment_0_column_type;
-	typedef common_column_type moment_2_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::moments_type > moments_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::moments_type > wind_moments_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::desi_type >    desi_column_type;
 private:
-	free_seeing_column_type       m_free_seeing_column;
-	isoplanatic_angle_column_type m_isoplanatic_angle_column;
-	altitude_eff_column_type      m_altitude_eff_column;
-	desi_time_constant_column_type m_desi_time_constant_column;
-	time_constant_column_type     m_time_constant_column;
-	moment_0_column_type          m_moment_0_column;
-	moment_2_column_type          m_moment_2_column;
+	moments_column_type m_moments_column;
+	wind_moments_column_type m_wind_moments_column;
+	desi_column_type m_desi_column;
 public:
 	storage( size_type reserve = 60 );
 
-	const free_seeing_column_type&       free_seeing()       const { return m_free_seeing_column; }
-	const isoplanatic_angle_column_type& isoplanatic_angle() const { return m_isoplanatic_angle_column; }
-	const altitude_eff_column_type&      altitude_eff()      const { return m_altitude_eff_column; }
-	const desi_time_constant_column_type& desi_time_constant() const { return m_desi_time_constant_column; }
-	const time_constant_column_type&     time_constant()     const { return m_time_constant_column; }
-	const moment_0_column_type&          moment_0()          const { return m_moment_0_column; }
-	const moment_2_column_type&          moment_2()          const { return m_moment_2_column; }
+	const moments_column_type& moments() const { return m_moments_column; }
+	const wind_moments_column_type& wind_moments() const { return m_wind_moments_column; }
+	const desi_column_type& desi() const { return m_desi_column; }
 
-	void push_back( const detail::sample& integrals ) {
-		table_storage< detail::sample >::push_back( integrals );
+	void push_back( const moments_type& moments, const moments_type& wind_moments, const desi_type& desi ) {
+		table_storage< detail::sample >::push_back( detail::sample( moments, wind_moments, desi ) );
 	}
 };
 
-namespace filter {
-	class regular {
-	public:
-		typedef atmos::integrals::detail::sample sample_type;
-	public:
-		bool operator() ( const sample_type& x ) const;
-	private:
-		bool is_regular( double x ) const;
-	};
-}
-
 class functor_factory;
 class functor;
 
@@ -121,7 +111,7 @@
 	private:
 		storage_type m_storage;
 
-		dcf_matrix_type calculate( const weight_type& wf, unsigned int npower, const double* powers );
+		dcf_matrix_type calculate( const weight_type& wf, const vector_type& powers );
 	public:
 		functor_factory();
 		functor_factory( const functor_factory& factory );
@@ -138,15 +128,18 @@
 
 		typedef atmos::indices::indices_type indices_type;
 		typedef atmos::indices::indesis_type indesis_type;
-		typedef atmos::integrals::detail::sample result_type;
+		typedef boost::tuple< atmos::integrals::moments_type, atmos::integrals::moments_type, atmos::integrals::desi_type > result_type;
 	private:
 		dcf_matrix_type m_indices_dcf;
 		dcf_matrix_type m_wind_indices_dcf;
+
+		const vector_type& m_powers;
+		const vector_type& m_wind_powers;
 	public:
-		functor( const dcf_matrix_type& indices_dcf, const dcf_matrix_type& wind_indices_dcf ):
-			m_indices_dcf( indices_dcf ), m_wind_indices_dcf( wind_indices_dcf ) {}
+		functor( const dcf_matrix_type& indices_dcf, const vector_type& powers, const dcf_matrix_type& wind_indices_dcf, const vector_type& wind_powers ):
+			m_indices_dcf( indices_dcf ), m_powers(powers), m_wind_indices_dcf( wind_indices_dcf ), m_wind_powers(wind_powers) {}
 
-		result_type operator() ( const indices_type& indices, const indesis_type& indesis, const indices_type& wind_indices, double mz ) const;
+		result_type operator() ( const indices_type& indices, const indices_type& wind_indices, const indesis_type& desi, double mz ) const;
 };
 
 }}

atmos-2.97.2.tar.gz/src/worksheet.cpp -> atmos-2.97.3.tar.gz/src/worksheet.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: worksheet.cpp 203 2011-02-03 17:40:12Z matwey $
+    $Id: worksheet.cpp 210 2011-03-02 18:59:06Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -219,9 +219,9 @@
 				double Jtrans_err = dimm::turbulence( context().m_pproxy.dimm_scale(), context().m_pproxy.dimm_Ktransversal(),  yerrs(1) );
 
 				dimm_turbulence = 0.5*( Jlong + Jtrans ) / m_mean_mz;
-				dimm_turbulence_err = 0.5 / m_mean_mz * std::sqrt( Jlong_err*Jlong_err + Jtrans_err*Jtrans_err );
+				dimm_turbulence_err = std::sqrt( Jlong_err*Jlong_err + Jtrans_err*Jtrans_err ) / (Jlong + Jtrans);
 				dimm_seeing = dimm::seeing( dimm_turbulence );
-				dimm_seeing_err = 0.6 * dimm_seeing / dimm_turbulence * dimm_turbulence_err;
+				dimm_seeing_err = 0.6*dimm_turbulence_err;
 			}
 
 			indices::flux_type flux_avg = average( m_indices.flux().begin(), m_indices.flux().end(), mean< indices::flux_type >() );
@@ -238,36 +238,39 @@
 
 			integrals::storage integrals;
 			for( indices::storage::const_iterator it = m_indices.begin(); it != m_indices.end(); ++it ) {
-				integrals.push_back( ifun( it->indices, it->indesis, it->wind_indices, m_mean_mz ) );
-			}
+				const boost::tuple<integrals::moments_type,integrals::moments_type,integrals::desi_type>& moments = ifun( it->indices, it->wind_indices, it->indesis, m_mean_mz );
 
-			integrals.unstable_remove_if( integrals::filter::regular() );
+				integrals.push_back( moments.get<0>(), moments.get<1>(), moments.get<2>() );
+			}
 
-			double free_seeing_avg       = average( integrals.free_seeing().begin(),       integrals.free_seeing().end(),       mean< double >() );
-			double isoplanatic_angle_avg = average( integrals.isoplanatic_angle().begin(), integrals.isoplanatic_angle().end(), mean< double >() );
-			double altitude_eff_avg      = average( integrals.altitude_eff().begin(),      integrals.altitude_eff().end(),      mean< double >() );
-			double desi_time_constant_avg= average( integrals.desi_time_constant().begin(),integrals.desi_time_constant().end(),mean< double >() );
-			double time_constant_avg     = average( integrals.time_constant().begin(),     integrals.time_constant().end(),     mean< double >() );
-			double moment_0_avg          = average( integrals.moment_0().begin(),          integrals.moment_0().end(),          mean< double >() );
-			double moment_2_avg          = average( integrals.moment_2().begin(),          integrals.moment_2().end(),          mean< double >() );
-			double free_seeing_err       = average( integrals.free_seeing().begin(),       integrals.free_seeing().end(),       error_of_mean< double >(free_seeing_avg) );
-			double isoplanatic_angle_err = average( integrals.isoplanatic_angle().begin(), integrals.isoplanatic_angle().end(), error_of_mean< double >(isoplanatic_angle_avg) );
-			double altitude_eff_err      = average( integrals.altitude_eff().begin(),      integrals.altitude_eff().end(),      error_of_mean< double >(altitude_eff_avg) );
-			double desi_time_constant_err= average( integrals.desi_time_constant().begin(),integrals.desi_time_constant().end(),error_of_mean< double >(desi_time_constant_avg) );
-			double time_constant_err     = average( integrals.time_constant().begin(),     integrals.time_constant().end(),     error_of_mean< double >(time_constant_avg) );
-			double moment_0_err          = average( integrals.moment_0().begin(),          integrals.moment_0().end(),          error_of_mean< double >(moment_0_avg) );
-			double moment_2_err          = average( integrals.moment_2().begin(),          integrals.moment_2().end(),          error_of_mean< double >(moment_2_avg) );
+			integrals::moments_type moments = average( integrals.moments().begin(), integrals.moments().end(), mean< integrals::moments_type >() );
+			integrals::moments_type moments_err = average( integrals.moments().begin(), integrals.moments().end(), error_of_mean< integrals::moments_type >(moments) );
+			integrals::moments_type wind_moments = average( integrals.wind_moments().begin(), integrals.wind_moments().end(), mean< integrals::moments_type >() );
+			integrals::moments_type wind_moments_err = average( integrals.wind_moments().begin(), integrals.wind_moments().end(), error_of_mean< integrals::moments_type >(wind_moments) );
+			integrals::desi_type desi = average( integrals.desi().begin(), integrals.desi().end(), mean< integrals::desi_type >() );
+			integrals::desi_type desi_err = average( integrals.desi().begin(), integrals.desi().end(), error_of_mean< integrals::desi_type >(desi) );
+
+			double moment_0 = moments(integrals::moment_0);
+			double moment_2 = moments(integrals::moment_2);
+			double moment_0_err = moments_err(integrals::moment_0)/moments(integrals::moment_0);
+			double moment_2_err = moments_err(integrals::moment_2)/moments(integrals::moment_2);
+			
+			std::pair<double,double> free_seeing = integrals::seeing(moment_0,moment_0_err);
+			std::pair<double,double> isoplanatic_angle = integrals::isoplanatic_angle(moments(integrals::moment_53),moments_err(integrals::moment_53)/moments(integrals::moment_53));
+			std::pair<double,double> altitude_eff = integrals::effective_altitude(moment_0,moment_0_err,moments(integrals::moment_1),moments_err(integrals::moment_1)/moments(integrals::moment_1));
+			std::pair<double,double> desi_time_constant = integrals::time_constant_desi(desi,desi_err/desi);
+			std::pair<double,double> time_constant = integrals::time_constant(wind_moments(integrals::moment_0),wind_moments_err(integrals::moment_0)/wind_moments(integrals::moment_0),moment_0,moment_0_err);
 
 			context().m_file_output.atmosphere( meantime,
-				free_seeing_avg, free_seeing_err,
+				free_seeing.first, free_seeing.second,
 				dimm_seeing, dimm_seeing_err,
 				dimm_turbulence, dimm_turbulence_err,
-				isoplanatic_angle_avg, isoplanatic_angle_err,
-				altitude_eff_avg, altitude_eff_err,
-				desi_time_constant_avg, desi_time_constant_err,
-				time_constant_avg, time_constant_err,
-				moment_0_avg, moment_0_err,
-				moment_2_avg, moment_2_err );
+				isoplanatic_angle.first, isoplanatic_angle.second,
+				altitude_eff.first, altitude_eff.second,
+				desi_time_constant.first, desi_time_constant.second,
+				time_constant.first, time_constant.second,
+				moment_0, moment_0_err,
+				moment_2, moment_2_err );
 
 			/* mixture */
 			mixture::storage mixture_samples( m_indices, m_dimm, m_dimm.size() > 1 ? context().m_pproxy.dimm_scale() : 0 );