Open Build Service

We truncated the diff of some files because they were too big. If you want to see the full diff for every file, click here.

Changes of Revision 9

atmos.spec Changed

atmos-2.96.9.tar.gz/configure -> atmos-2.97.1.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.96.9.
+# Generated by GNU Autoconf 2.63 for atmos 2.97.1.
 #
 # 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.96.9'
-PACKAGE_STRING='atmos 2.96.9'
+PACKAGE_VERSION='2.97.1'
+PACKAGE_STRING='atmos 2.97.1'
 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.96.9 to adapt to many kinds of systems.
+\`configure' configures atmos 2.97.1 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.96.9:";;
+     short | recursive ) echo "Configuration of atmos 2.97.1:";;
    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.96.9
+atmos configure 2.97.1
 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.96.9, which was
+It was created by atmos $as_me 2.97.1, 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.96.9'
+ VERSION='2.97.1'
 
 
 cat >>confdefs.h <<_ACEOF
@@ -3978,11 +3978,15 @@
 
 
 
-for ac_header in boost/numeric/ublas/io.hpp \
+
+
+for ac_header in boost/numeric/ublas/banded.hpp \
+	boost/numeric/ublas/matrix_expression.hpp \
 	boost/numeric/ublas/matrix.hpp \
 	boost/numeric/ublas/matrix_proxy.hpp \
 	boost/numeric/ublas/symmetric.hpp \
 	boost/numeric/ublas/triangular.hpp \
+	boost/numeric/ublas/vector_expression.hpp \
 	boost/numeric/ublas/vector.hpp \
 	boost/numeric/ublas/vector_proxy.hpp
 do
@@ -4152,6 +4156,10 @@
 
 
 
+
+
+
+
 for ac_header in boost/algorithm/string.hpp \
 	boost/bind.hpp \
 	boost/cstdint.hpp \
@@ -4161,10 +4169,14 @@
 	boost/format.hpp \
 	boost/function.hpp \
 	boost/io/ios_state.hpp \
-	boost/iterator.hpp \
+	boost/iterator/transform_iterator.hpp \
 	boost/lexical_cast.hpp \
+	boost/math/tools/roots.hpp \
+	boost/mem_fn.hpp \
+	boost/mpl/if.hpp \
 	boost/program_options.hpp \
-	boost/tuple/tuple.hpp
+	boost/tuple/tuple.hpp \
+	boost/type_traits/is_complex.hpp
 do
 as_ac_Header=`$as_echo "ac_cv_header_$ac_header" | $as_tr_sh`
 if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
@@ -4321,160 +4333,6 @@
 
 
 
-for ac_header in boost/math/tools/roots.hpp \
-	tr1/tuple
-do
-as_ac_Header=`$as_echo "ac_cv_header_$ac_header" | $as_tr_sh`
-if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
-  { $as_echo "$as_me:$LINENO: checking for $ac_header" >&5
-$as_echo_n "checking for $ac_header... " >&6; }
-if { as_var=$as_ac_Header; eval "test \"\${$as_var+set}\" = set"; }; then
-  $as_echo_n "(cached) " >&6
-fi
-ac_res=`eval 'as_val=${'$as_ac_Header'}
-		 $as_echo "$as_val"'`
-	       { $as_echo "$as_me:$LINENO: result: $ac_res" >&5
-$as_echo "$ac_res" >&6; }
-else
-  # Is the header compilable?
-{ $as_echo "$as_me:$LINENO: checking $ac_header usability" >&5
-$as_echo_n "checking $ac_header usability... " >&6; }
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-$ac_includes_default
-#include <$ac_header>
-_ACEOF
-rm -f conftest.$ac_objext
-if { (ac_try="$ac_compile"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval ac_try_echo="\"\$as_me:$LINENO: $ac_try_echo\""
-$as_echo "$ac_try_echo") >&5
-  (eval "$ac_compile") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  $as_echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } && {
-	 test -z "$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       } && test -s conftest.$ac_objext; then
-  ac_header_compiler=yes
-else
-  $as_echo "$as_me: failed program was:" >&5
-sed 's/^/| /' conftest.$ac_ext >&5
-
-	ac_header_compiler=no
-fi
-
-rm -f core conftest.err conftest.$ac_objext conftest.$ac_ext
-{ $as_echo "$as_me:$LINENO: result: $ac_header_compiler" >&5
-$as_echo "$ac_header_compiler" >&6; }
-
-# Is the header present?
-{ $as_echo "$as_me:$LINENO: checking $ac_header presence" >&5
-$as_echo_n "checking $ac_header presence... " >&6; }
-cat >conftest.$ac_ext <<_ACEOF
-/* confdefs.h.  */
-_ACEOF
-cat confdefs.h >>conftest.$ac_ext
-cat >>conftest.$ac_ext <<_ACEOF
-/* end confdefs.h.  */
-#include <$ac_header>
-_ACEOF
-if { (ac_try="$ac_cpp conftest.$ac_ext"
-case "(($ac_try" in
-  *\"* | *\`* | *\\*) ac_try_echo=\$ac_try;;
-  *) ac_try_echo=$ac_try;;
-esac
-eval ac_try_echo="\"\$as_me:$LINENO: $ac_try_echo\""
-$as_echo "$ac_try_echo") >&5
-  (eval "$ac_cpp conftest.$ac_ext") 2>conftest.er1
-  ac_status=$?
-  grep -v '^ *+' conftest.er1 >conftest.err
-  rm -f conftest.er1
-  cat conftest.err >&5
-  $as_echo "$as_me:$LINENO: \$? = $ac_status" >&5
-  (exit $ac_status); } >/dev/null && {
-	 test -z "$ac_cxx_preproc_warn_flag$ac_cxx_werror_flag" ||
-	 test ! -s conftest.err
-       }; then
-  ac_header_preproc=yes
-else

atmos-2.96.9.tar.gz/configure.in -> atmos-2.97.1.tar.gz/configure.in Changed

@@ -2,7 +2,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ(2.59)
-AC_INIT(atmos, 2.96.9, mass-general@curl.sai.msu.ru)
+AC_INIT(atmos, 2.97.1, mass-general@curl.sai.msu.ru)
 AM_INIT_AUTOMAKE()
 AC_CONFIG_HEADERS(src/config.h)
 AC_CONFIG_SRCDIR(src/atmos.h)
@@ -18,11 +18,13 @@
 
 # Checks for header files.
 AC_CHECK_HEADERS(stdlib.h string.h unistd.h)
-AC_CHECK_HEADERS(boost/numeric/ublas/io.hpp \
+AC_CHECK_HEADERS(boost/numeric/ublas/banded.hpp \
+	boost/numeric/ublas/matrix_expression.hpp \
 	boost/numeric/ublas/matrix.hpp \
 	boost/numeric/ublas/matrix_proxy.hpp \
 	boost/numeric/ublas/symmetric.hpp \
 	boost/numeric/ublas/triangular.hpp \
+	boost/numeric/ublas/vector_expression.hpp \
 	boost/numeric/ublas/vector.hpp \
 	boost/numeric/ublas/vector_proxy.hpp,,AC_MSG_ERROR(can't find boost ublas headers))
 
@@ -35,13 +37,14 @@
 	boost/format.hpp \
 	boost/function.hpp \
 	boost/io/ios_state.hpp \
-	boost/iterator.hpp \
+	boost/iterator/transform_iterator.hpp \
 	boost/lexical_cast.hpp \
+	boost/math/tools/roots.hpp \
+	boost/mem_fn.hpp \
+	boost/mpl/if.hpp \
 	boost/program_options.hpp \
-	boost/tuple/tuple.hpp,,AC_MSG_ERROR(can't find boost headers))
-
-AC_CHECK_HEADERS(boost/math/tools/roots.hpp \
-	tr1/tuple)
+	boost/tuple/tuple.hpp \
+	boost/type_traits/is_complex.hpp,,AC_MSG_ERROR(can't find boost headers))
 
 AC_CHECK_HEADERS(gsl/gsl_errno.h \
 	gsl/gsl_integration.h ,,AC_MSG_ERROR(can't find gsl headers))

atmos-2.96.9.tar.gz/src/Makefile.am -> atmos-2.97.1.tar.gz/src/Makefile.am Changed

@@ -2,15 +2,48 @@
 
 atmos_CPPFLAGS = -I @abs_top_srcdir@/third-party/lsp/include -I @abs_top_srcdir@/third-party/yalpa/include
 
-atmos_SOURCES  = apply_to_all.h atmos.h average.h background_factory.h background_factory.cpp \
-		dimm.h dimm.cpp \
-		file_input.cpp file_input.h file_output.cpp file_output.h fft.h \
-		indices.cpp indices.h integrals.cpp integrals.h \
-		io.cpp io.h gsl_wrap.h lambert.h main.cpp nocase_less.h params.h \
-		photometry.h profile.cpp profile.h spectral_factory.cpp \
-		spectral_factory.h spline.h table_storage.h typesdef.h unstable_remove_if.h utils.h \
-		vector_generator.h worksheet.h worksheet.cpp \
-		weif.cpp weif.h
+atmos_SOURCES  = apply_to_all.h \
+	atmos.h \
+	average.h \
+	background_factory.h \
+	covariance_decomposition.h \
+	dimm.h \
+	fft.h \
+	file_input.h \
+	file_output.h \
+	gsl_wrap.h \
+	indices.h \
+	integrals.h \
+	io.h \
+	lambert.h \
+	local_context.h \
+	nocase_less.h \
+	params.h \
+	photometry.h \
+	profile.h \
+	spectral_factory.h \
+	spline.h \
+	table_storage.h \
+	typesdef.h \
+	unstable_remove_if.h \
+	utils.h \
+	vector_generator.h \
+	weif.h \
+	worksheet.h \
+	background_factory.cpp \
+	dimm.cpp \
+	file_input.cpp \
+	file_output.cpp \
+	indices.cpp \
+	integrals.cpp \
+	io.cpp \
+	main.cpp \
+	params.cpp \
+	profile.cpp \
+	spectral_factory.cpp \
+	utils.cpp \
+	weif.cpp \
+	worksheet.cpp
 
 atmos_LDADD    = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
 atmos_DEPENDENCIES = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a

atmos-2.96.9.tar.gz/src/Makefile.in -> atmos-2.97.1.tar.gz/src/Makefile.in Changed

@@ -50,9 +50,10 @@
 	atmos-dimm.$(OBJEXT) atmos-file_input.$(OBJEXT) \
 	atmos-file_output.$(OBJEXT) atmos-indices.$(OBJEXT) \
 	atmos-integrals.$(OBJEXT) atmos-io.$(OBJEXT) \
-	atmos-main.$(OBJEXT) atmos-profile.$(OBJEXT) \
-	atmos-spectral_factory.$(OBJEXT) atmos-worksheet.$(OBJEXT) \
-	atmos-weif.$(OBJEXT)
+	atmos-main.$(OBJEXT) atmos-params.$(OBJEXT) \
+	atmos-profile.$(OBJEXT) atmos-spectral_factory.$(OBJEXT) \
+	atmos-utils.$(OBJEXT) atmos-weif.$(OBJEXT) \
+	atmos-worksheet.$(OBJEXT)
 atmos_OBJECTS = $(am_atmos_OBJECTS)
 DEFAULT_INCLUDES = -I.@am__isrc@
 depcomp = $(SHELL) $(top_srcdir)/depcomp
@@ -159,15 +160,48 @@
 top_builddir = @top_builddir@
 top_srcdir = @top_srcdir@
 atmos_CPPFLAGS = -I @abs_top_srcdir@/third-party/lsp/include -I @abs_top_srcdir@/third-party/yalpa/include
-atmos_SOURCES = apply_to_all.h atmos.h average.h background_factory.h background_factory.cpp \
-		dimm.h dimm.cpp \
-		file_input.cpp file_input.h file_output.cpp file_output.h fft.h \
-		indices.cpp indices.h integrals.cpp integrals.h \
-		io.cpp io.h gsl_wrap.h lambert.h main.cpp nocase_less.h params.h \
-		photometry.h profile.cpp profile.h spectral_factory.cpp \
-		spectral_factory.h spline.h table_storage.h typesdef.h unstable_remove_if.h utils.h \
-		vector_generator.h worksheet.h worksheet.cpp \
-		weif.cpp weif.h
+atmos_SOURCES = apply_to_all.h \
+	atmos.h \
+	average.h \
+	background_factory.h \
+	covariance_decomposition.h \
+	dimm.h \
+	fft.h \
+	file_input.h \
+	file_output.h \
+	gsl_wrap.h \
+	indices.h \
+	integrals.h \
+	io.h \
+	lambert.h \
+	local_context.h \
+	nocase_less.h \
+	params.h \
+	photometry.h \
+	profile.h \
+	spectral_factory.h \
+	spline.h \
+	table_storage.h \
+	typesdef.h \
+	unstable_remove_if.h \
+	utils.h \
+	vector_generator.h \
+	weif.h \
+	worksheet.h \
+	background_factory.cpp \
+	dimm.cpp \
+	file_input.cpp \
+	file_output.cpp \
+	indices.cpp \
+	integrals.cpp \
+	io.cpp \
+	main.cpp \
+	params.cpp \
+	profile.cpp \
+	spectral_factory.cpp \
+	utils.cpp \
+	weif.cpp \
+	worksheet.cpp
 
 atmos_LDADD = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
 atmos_DEPENDENCIES = @abs_top_builddir@/third-party/yalpa/src/.libs/libyalpa.a
@@ -278,8 +312,10 @@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-integrals.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-io.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-main.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-params.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-profile.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-spectral_factory.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-utils.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-weif.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/atmos-worksheet.Po@am__quote@
 
@@ -409,6 +445,20 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-main.obj `if test -f 'main.cpp'; then $(CYGPATH_W) 'main.cpp'; else $(CYGPATH_W) '$(srcdir)/main.cpp'; fi`
 
+atmos-params.o: params.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-params.o -MD -MP -MF $(DEPDIR)/atmos-params.Tpo -c -o atmos-params.o `test -f 'params.cpp' || echo '$(srcdir)/'`params.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-params.Tpo $(DEPDIR)/atmos-params.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='params.cpp' object='atmos-params.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-params.o `test -f 'params.cpp' || echo '$(srcdir)/'`params.cpp
+
+atmos-params.obj: params.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-params.obj -MD -MP -MF $(DEPDIR)/atmos-params.Tpo -c -o atmos-params.obj `if test -f 'params.cpp'; then $(CYGPATH_W) 'params.cpp'; else $(CYGPATH_W) '$(srcdir)/params.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-params.Tpo $(DEPDIR)/atmos-params.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='params.cpp' object='atmos-params.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-params.obj `if test -f 'params.cpp'; then $(CYGPATH_W) 'params.cpp'; else $(CYGPATH_W) '$(srcdir)/params.cpp'; fi`
+
 atmos-profile.o: profile.cpp
 @am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-profile.o -MD -MP -MF $(DEPDIR)/atmos-profile.Tpo -c -o atmos-profile.o `test -f 'profile.cpp' || echo '$(srcdir)/'`profile.cpp
 @am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-profile.Tpo $(DEPDIR)/atmos-profile.Po
@@ -437,19 +487,19 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-spectral_factory.obj `if test -f 'spectral_factory.cpp'; then $(CYGPATH_W) 'spectral_factory.cpp'; else $(CYGPATH_W) '$(srcdir)/spectral_factory.cpp'; fi`
 
-atmos-worksheet.o: worksheet.cpp
-@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.o -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
-@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.o' libtool=no @AMDEPBACKSLASH@
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
-
-atmos-worksheet.obj: worksheet.cpp
-@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.obj -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
-@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
-@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.obj' libtool=no @AMDEPBACKSLASH@
+atmos-utils.o: utils.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-utils.o -MD -MP -MF $(DEPDIR)/atmos-utils.Tpo -c -o atmos-utils.o `test -f 'utils.cpp' || echo '$(srcdir)/'`utils.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-utils.Tpo $(DEPDIR)/atmos-utils.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='utils.cpp' object='atmos-utils.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-utils.o `test -f 'utils.cpp' || echo '$(srcdir)/'`utils.cpp
+
+atmos-utils.obj: utils.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-utils.obj -MD -MP -MF $(DEPDIR)/atmos-utils.Tpo -c -o atmos-utils.obj `if test -f 'utils.cpp'; then $(CYGPATH_W) 'utils.cpp'; else $(CYGPATH_W) '$(srcdir)/utils.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-utils.Tpo $(DEPDIR)/atmos-utils.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='utils.cpp' object='atmos-utils.obj' libtool=no @AMDEPBACKSLASH@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
-@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-utils.obj `if test -f 'utils.cpp'; then $(CYGPATH_W) 'utils.cpp'; else $(CYGPATH_W) '$(srcdir)/utils.cpp'; fi`
 
 atmos-weif.o: weif.cpp
 @am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-weif.o -MD -MP -MF $(DEPDIR)/atmos-weif.Tpo -c -o atmos-weif.o `test -f 'weif.cpp' || echo '$(srcdir)/'`weif.cpp
@@ -465,6 +515,20 @@
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-weif.obj `if test -f 'weif.cpp'; then $(CYGPATH_W) 'weif.cpp'; else $(CYGPATH_W) '$(srcdir)/weif.cpp'; fi`
 
+atmos-worksheet.o: worksheet.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.o -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.o' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.o `test -f 'worksheet.cpp' || echo '$(srcdir)/'`worksheet.cpp
+
+atmos-worksheet.obj: worksheet.cpp
+@am__fastdepCXX_TRUE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -MT atmos-worksheet.obj -MD -MP -MF $(DEPDIR)/atmos-worksheet.Tpo -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+@am__fastdepCXX_TRUE@	$(am__mv) $(DEPDIR)/atmos-worksheet.Tpo $(DEPDIR)/atmos-worksheet.Po
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	source='worksheet.cpp' object='atmos-worksheet.obj' libtool=no @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(atmos_CPPFLAGS) $(CPPFLAGS) $(AM_CXXFLAGS) $(CXXFLAGS) -c -o atmos-worksheet.obj `if test -f 'worksheet.cpp'; then $(CYGPATH_W) 'worksheet.cpp'; else $(CYGPATH_W) '$(srcdir)/worksheet.cpp'; fi`
+
 ID: $(HEADERS) $(SOURCES) $(LISP) $(TAGS_FILES)
 	list='$(SOURCES) $(HEADERS) $(LISP) $(TAGS_FILES)'; \
 	unique=`for i in $$list; do \

atmos-2.96.9.tar.gz/src/apply_to_all.h -> atmos-2.97.1.tar.gz/src/apply_to_all.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: apply_to_all.h 155 2010-09-19 08:30:20Z matwey $
+    $Id: apply_to_all.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -21,10 +21,10 @@
 
 #include <cmath>
 
-#include "lambert.h"
+#include <boost/numeric/ublas/matrix_expression.hpp>
+#include <boost/numeric/ublas/vector_expression.hpp>
 
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/vector.hpp>
+#include "lambert.h"
 
 namespace utils {
 namespace ublas {
@@ -78,6 +78,14 @@
 			return x + value_type(1);
 		}
 	};
+	template <class T> class abs {
+	public:
+		typedef T value_type;
+		typedef T result_type;
+		static result_type apply( const value_type& x) {
+			return std::abs(x);
+		}
+	};
 }
 
 using namespace boost::numeric::ublas;
@@ -146,6 +154,15 @@
 	return apply_to_all( e, detail::lambert_W0< value_type >() );
 }
 
+template<class E> BOOST_UBLAS_INLINE typename vector_unary_traits<E, detail::abs< typename E::value_type > >::result_type abs(const vector_expression<E> &e ) {
+	typedef typename E::value_type value_type;
+	return apply_to_all( e, detail::abs< value_type >() );
+}
+template<class E> BOOST_UBLAS_INLINE typename matrix_unary1_traits<E, detail::abs< typename E::value_type > >::result_type abs(const matrix_expression<E> &e ) {
+	typedef typename E::value_type value_type;
+	return apply_to_all( e, detail::abs< value_type >() );
+}
+
 } //ublas
 } //utils

atmos-2.96.9.tar.gz/src/atmos.h -> atmos-2.97.1.tar.gz/src/atmos.h Changed

atmos-2.96.9.tar.gz/src/average.h -> atmos-2.97.1.tar.gz/src/average.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: average.h 143 2010-04-04 08:42:34Z matwey $
+    $Id: average.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -20,153 +20,156 @@
 #define _AVERAGE_H
 
 #include <cmath>
-#include <iterator>
+#include <iterator> // std::iterator_traits
 
-#include <boost/numeric/ublas/matrix.hpp>
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/io.hpp>
+#include <boost/numeric/ublas/matrix_expression.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/vector_expression.hpp>
 
 #include "apply_to_all.h"
 
 namespace utils {
 namespace algo {
 
-	template<class T> struct mean;
-	template<class T> struct variance;
+namespace detail {
 
-	template<class T> struct mean {
-	public:
-		typedef T value_type;
-		typedef T result_type;
-
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += x;
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( x );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size;
-		}
-	private:
-	};
-
-	struct variance_traits {
-		static double prod( double x, double y ){
-			return x*y;
-		}
-		static float prod( float x, float y ){
-			return x*y;
-		}
-		template< class V1, class V2 > static
-		typename boost::numeric::ublas::vector_binary_traits< V1, V2, boost::numeric::ublas::scalar_multiplies< typename V1::value_type, typename V2::value_type > >::result_type
-		prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
-			return element_prod( e1, e2 );
-		}
-		template< class M1, class M2 > static
-		typename boost::numeric::ublas::matrix_binary_traits< M1, M2, boost::numeric::ublas::scalar_multiplies< typename M1::value_type, typename M2::value_type > >::result_type
-		prod( const boost::numeric::ublas::matrix_expression< M1 >& e1, const boost::numeric::ublas::matrix_expression< M2 >& e2 ) {
-			return element_prod( e1, e2 );
-		}
-
-		static double sqrt( double x ){
-			return std::sqrt(x);
-		}
-		static float sqrt( float x ){
-			return std::sqrt(x);
-		}
-		template< class V1 > static
-		typename boost::numeric::ublas::vector_unary_traits< V1, utils::ublas::detail::sqrt< typename V1::value_type > >::result_type
-		sqrt( const boost::numeric::ublas::vector_expression< V1 >& e1) {
-			return utils::ublas::sqrt( e1 );
-		}
-		template< class M1 > static
-		typename boost::numeric::ublas::matrix_unary1_traits< M1, utils::ublas::detail::sqrt< typename M1::value_type > >::result_type
-		sqrt( const boost::numeric::ublas::matrix_expression< M1 >& e1) {
-			return utils::ublas::sqrt( e1 );
-		}
-	};
-
-	struct covariance_traits {
-		template< class V1, class V2 > static
-		boost::numeric::ublas::symmetric_adaptor< const typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type, boost::numeric::ublas::upper >
-		prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
-			typedef typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type outer_prod_type;
-			return boost::numeric::ublas::symmetric_adaptor< const outer_prod_type, boost::numeric::ublas::upper >( outer_prod( e1, e2 ) );
-		}
-	};
-	
-	template<class T> struct variance {
-	public:
-		typedef T value_type;
-		typedef T result_type;
-		
-		variance( const value_type& mean ): mean_(mean) {}
-		
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += variance_traits::prod( x - mean_, x - mean_ );
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( variance_traits::prod( x - mean_, x - mean_ ) );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size-1;
-		}
-	private:
-		value_type mean_;
-	};
-
-	template<class T> struct covariance {
-	public:
-		typedef T value_type;
-		typedef boost::numeric::ublas::symmetric_matrix< typename value_type::value_type > result_type;
-
-		covariance( const value_type& mean ): mean_(mean) {}
-
-		void operator() ( result_type& result, const value_type& x ) const {
-			result += covariance_traits::prod( x - mean_, x - mean_ );
-		}
-		result_type initial( const value_type& x ) const {
-			return result_type( covariance_traits::prod( x - mean_, x - mean_ ) );
-		}
-		template<class S> void result( result_type& result, S size ) const {
-			result /= size-1;
-		}
-	private:
-		value_type mean_;
-	};
-
-	template<class T> struct error_of_mean: variance<T> {
-		error_of_mean( const typename variance<T>::value_type& mean ): variance<T>(mean) {}
-		template<class S> void result( typename variance<T>::result_type& result, S size ) const {
-			variance<T>::result( result, size );
-			result = variance_traits::sqrt( result / size );
-		}
-	};
+struct variance_traits {
+	static double prod( double x, double y ){
+		return x*y;
+	}
+	static float prod( float x, float y ){
+		return x*y;
+	}
+	template< class V1, class V2 > static
+	typename boost::numeric::ublas::vector_binary_traits< V1, V2, boost::numeric::ublas::scalar_multiplies< typename V1::value_type, typename V2::value_type > >::result_type
+	prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
+		return element_prod( e1, e2 );
+	}
+	template< class M1, class M2 > static
+	typename boost::numeric::ublas::matrix_binary_traits< M1, M2, boost::numeric::ublas::scalar_multiplies< typename M1::value_type, typename M2::value_type > >::result_type
+	prod( const boost::numeric::ublas::matrix_expression< M1 >& e1, const boost::numeric::ublas::matrix_expression< M2 >& e2 ) {
+		return element_prod( e1, e2 );
+	}
+};
+
+struct error_of_mean_traits {
+	static double sqrt( double x ){
+		return std::sqrt(x);
+	}
+	static float sqrt( float x ){
+		return std::sqrt(x);
+	}
+	template< class V1 > static
+	typename boost::numeric::ublas::vector_unary_traits< V1, utils::ublas::detail::sqrt< typename V1::value_type > >::result_type
+	sqrt( const boost::numeric::ublas::vector_expression< V1 >& e1) {
+		return utils::ublas::sqrt( e1 );
+	}
+	template< class M1 > static
+	typename boost::numeric::ublas::matrix_unary1_traits< M1, utils::ublas::detail::sqrt< typename M1::value_type > >::result_type
+	sqrt( const boost::numeric::ublas::matrix_expression< M1 >& e1) {
+		return utils::ublas::sqrt( e1 );
+	}
+};
+
+struct covariance_traits {
+	template< class V1, class V2 > static
+	boost::numeric::ublas::symmetric_adaptor< const typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type, boost::numeric::ublas::upper >
+	prod( const boost::numeric::ublas::vector_expression< V1 >& e1, const boost::numeric::ublas::vector_expression< V2 >& e2 ) {
+		typedef typename boost::numeric::ublas::vector_matrix_binary_traits<V1, V2, boost::numeric::ublas::scalar_multiplies<typename V1::value_type, typename V2::value_type> >::result_type outer_prod_type;
+		return boost::numeric::ublas::symmetric_adaptor< const outer_prod_type, boost::numeric::ublas::upper >( outer_prod( e1, e2 ) );
+	}
+};
 	
-	template< class Iterator, template<class T> class F >
-	typename F< typename std::iterator_traits<Iterator>::value_type >::result_type
-	average( Iterator begin, Iterator end, F< typename std::iterator_traits<Iterator>::value_type > op ) {
-		typedef Iterator iterator;
-		typedef typename std::iterator_traits< iterator >::value_type value_type;
-		typedef typename std::iterator_traits< iterator >::difference_type difference_type;
-		typedef difference_type size_type;
-		typedef F< value_type > operator_type;
-		typedef typename operator_type::result_type result_type;
+}
+
+template<class T> struct mean {
+public:
+	typedef T value_type;

atmos-2.96.9.tar.gz/src/background_factory.cpp -> atmos-2.97.1.tar.gz/src/background_factory.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: background_factory.cpp 147 2010-04-17 09:25:25Z matwey $
+    $Id: background_factory.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,73 +16,79 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <stdexcept>
+#include <boost/format.hpp>
+
 #include "background_factory.h"
 
-#include <boost/format.hpp>
+namespace atmos {
+namespace background {
 
-background_factory::background_factory(const select_model_function_type& select_model):
+factory::factory(const select_model_function_type& select_model):
 	m_curmodel(0),
 	m_select_model(select_model) {
 }
-background_factory::background_factory(const background_factory& factory) {
+factory::factory(const factory& factory) {
 	m_select_model = factory.m_select_model;
 	m_curmodel = std::auto_ptr<model_type>(factory.m_curmodel->clone());
 }
-void background_factory::swap(background_factory& factory) throw() {
-	factory.m_select_model.swap( m_select_model );
-	std::swap( factory.m_curmodel, m_curmodel );
+void factory::swap(factory& f) throw() {
+	f.m_select_model.swap( m_select_model );
+	std::swap( f.m_curmodel, m_curmodel );
 }
-background_factory& background_factory::operator=(const background_factory& factory){
-	background_factory tmp( factory );
+factory& factory::operator=(const factory& f){
+	factory tmp( f );
 	tmp.swap(*this );
 	return *this;
 }
-background_factory::~background_factory() {
+factory::~factory() {
 }
-background_factory::vector_type background_factory::background( const boost::posix_time::ptime& time ) const {
+factory::flux_type factory::background( const time_type& time ) const {
 	if( ! m_curmodel.get() )
 		throw std::logic_error("There is not any background model");
 	return m_curmodel->background( time );
 }
-background_factory::model_type* background_factory::create_new_model( const boost::posix_time::ptime& time ) const {
+factory::model_type* factory::create_new_model( const time_type& time ) const {
 	return m_select_model( time );
 }
-background_factory::model_type* background_factory::release_current_model() {
+factory::model_type* factory::release_current_model() {
 	return m_curmodel.release();
 }
-void background_factory::handle_new_model( model_type* model ) {
+void factory::handle_new_model( model_type* model ) {
 	m_curmodel = std::auto_ptr<model_type>( model );
 }
 
-background_model::background_model( const boost::posix_time::ptime& time ): m_inittime(time) {
+namespace model {
+
+abstract::abstract( const time_type& time ): m_inittime(time) {
 }
-background_model::~background_model() {
+abstract::~abstract() {
 }
 
-background_daytime_model::background_daytime_model( const boost::posix_time::ptime& time ):
-	background_model(time),
+daytime::daytime( const time_type& time ):
+	abstract(time),
 	m_size(0) {
 }
-background_daytime_model::size_type background_daytime_model::size() const {
+daytime::size_type daytime::size() const {
 	return m_size;
 }
-background_daytime_model::vector_type background_daytime_model::background( const boost::posix_time::ptime& time ) const {
+daytime::flux_type daytime::background( const time_type& time ) const {
 	throw std::logic_error("Daytime background observation");
 }
-void background_daytime_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void daytime::accumulate( const time_type& time, const flux_type& flux ) {
 	m_size++;
 }
-background_model* background_daytime_model::clone() const {
-	return new background_daytime_model(*this);
+abstract* daytime::clone() const {
+	return new daytime(*this);
 }
 
-background_twilight_model::background_twilight_model( const boost::posix_time::ptime& time, const sun_altitude_function_type& sun_altitude ):
-	background_model(time),
+twilight::twilight( const time_type& time, const sun_altitude_function_type& sun_altitude ):
+	abstract(time),
 	m_size(0),
 	m_background(0),
 	m_sun_altitude(sun_altitude) {
 }
-background_twilight_model::value_type background_twilight_model::extrapolate( const boost::posix_time::ptime& timef ) const {
+value_type twilight::extrapolate( const time_type& timef ) const {
 	value_type hi = m_sun_altitude( inittime() ) * 180.0 / M_PI; // degrees
 	value_type hf = m_sun_altitude( timef ) * 180.0 / M_PI; // degrees
 	if( hi < -12.0 ) hi = -12.0;
@@ -90,62 +96,66 @@
 	value_type scale = 0.52189 - 0.26242*hi - 0.0235464*hi*hi; // model
 	return std::exp( scale * ( hf - hi ) );
 }
-background_twilight_model::size_type background_twilight_model::size() const {
+twilight::size_type twilight::size() const {
 	return m_size;
 }
-background_twilight_model::vector_type background_twilight_model::background( const boost::posix_time::ptime& time ) const {
+twilight::flux_type twilight::background( const time_type& time ) const {
 	return m_background * extrapolate( time );
 }
-void background_twilight_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void twilight::accumulate( const time_type& time, const flux_type& flux ) {
 
 	if( m_size == 0 ) {
-		m_background = mean / extrapolate( time );
+		m_background = flux / extrapolate( time );
 		m_size++;
 		return;
 	}
 
-	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + mean / extrapolate( time ) * ( value_type(1) / value_type( m_size + 1 ) );
+	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + flux / extrapolate( time ) * ( value_type(1) / value_type( m_size + 1 ) );
 	m_size++;
 }
-background_model* background_twilight_model::clone() const {
-	return new background_twilight_model(*this);
+abstract* twilight::clone() const {
+	return new twilight(*this);
 }
 
-background_night_model::background_night_model( const boost::posix_time::ptime& time ):
-	background_model(time),
+night::night( const time_type& time ):
+	abstract(time),
 	m_size(0),
 	m_background(0) {
 }
-background_night_model::size_type background_night_model::size() const {
+night::size_type night::size() const {
 	return m_size;
 }
-background_night_model::vector_type background_night_model::background( const boost::posix_time::ptime& time ) const {
+night::flux_type night::background( const time_type& time ) const {
 	return m_background;
 }
-void background_night_model::accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) {
+void night::accumulate( const time_type& time, const flux_type& flux ) {
 	if( m_size == 0 ) {
-		m_background = mean;
+		m_background = flux;
 		m_size++;
 		return;
 	}
 
-	assert( mean.size() == m_background.size() );
+	assert( flux.size() == m_background.size() );
 
-	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + mean * ( value_type(1) / value_type( m_size + 1 ) );
+	m_background = m_background * ( value_type( m_size ) / value_type( m_size + 1 ) ) + flux * ( value_type(1) / value_type( m_size + 1 ) );
 	m_size++;
 }
-background_model* background_night_model::clone() const {
-	return new background_night_model(*this);
+abstract* night::clone() const {
+	return new night(*this);
 }
 
-background_nonpoisson_filter::background_nonpoisson_filter( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp ):
+}
+
+namespace filter {
+
+nonpoisson::nonpoisson( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp ):
 	m_threshold( threshold ),
 	m_nonpoissonD( nonpoissonD ),
 	m_numexp( numexp ) {
 }
-background_nonpoisson_filter::~background_nonpoisson_filter() {
+nonpoisson::~nonpoisson() {
 }
-bool background_nonpoisson_filter::operator() ( const boost::posix_time::ptime& time, const flux_type& flux, const vars_type& vars ) const {
+bool nonpoisson::operator() ( const time_type& time, const flux_type& flux, const vars_type& vars ) const {
 	double poisson = ( vars(3,3)*flux(3) );
 	double relepsp = 2.0 * ( 1.0 + 1.0/flux(3) ) / ( m_numexp );
 	double pmax = m_nonpoissonD * ( 1.0 + m_threshold * std::sqrt(relepsp) );
@@ -155,3 +165,6 @@

atmos-2.96.9.tar.gz/src/background_factory.h -> atmos-2.97.1.tar.gz/src/background_factory.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: background_factory.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: background_factory.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,8 +19,7 @@
 #ifndef _BACKGROUND_FACTORY_H
 #define _BACKGROUND_FACTORY_H
 
-#include <list>
-#include <memory>
+#include <memory> // std::auto_ptr
 
 #include <boost/function.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
@@ -28,133 +27,137 @@
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 
-using namespace boost::numeric;
+namespace atmos {
+namespace background {
 
-/* Forward declaration */
-class background_factory;
+using namespace boost::numeric::ublas;
 
-class background_model;
-class background_daytime_model;
-class background_twilight_model;
-class background_night_model;
-
-class background_abstract_filter;
-class background_nonpoisson_filter;
-
-/* Declaration */
-class background_factory {
-	public:
-		typedef background_abstract_filter filter_type;
-		typedef background_model  model_type;
-
-		typedef boost::function1<model_type*, const boost::posix_time::ptime&> select_model_function_type;
-
-		typedef double                     value_type;
-		typedef ublas::vector<value_type>  vector_type;
-		typedef vector_type::size_type     size_type;
-
-		typedef vector_type flux_type;
-		typedef ublas::symmetric_matrix<value_type> vars_type;
-	private:
-		std::auto_ptr< model_type > m_curmodel;
-		select_model_function_type  m_select_model;
-	public:
-		background_factory(const select_model_function_type& select_model);
-		background_factory(const background_factory& factory);
-		void swap(background_factory& factory) throw();
-		background_factory& operator=(const background_factory& factory);
-		~background_factory();
-
-		vector_type background( const boost::posix_time::ptime& time ) const;
-
-		model_type* create_new_model( const boost::posix_time::ptime& time ) const;
-		model_type* release_current_model();
-		void handle_new_model( model_type* model );
-};
-
-class background_model {
-	public:
-		typedef background_factory::value_type  value_type;
-		typedef background_factory::vector_type vector_type;
-		typedef background_factory::size_type   size_type;
-	private:
-		boost::posix_time::ptime m_inittime;
-	public:
-		background_model( const boost::posix_time::ptime& time );
-		virtual ~background_model() = 0;
-
-		const boost::posix_time::ptime& inittime() const { return m_inittime; }
-
-		virtual size_type size() const = 0;
-		virtual vector_type background( const boost::posix_time::ptime& time ) const = 0;
-		virtual void accumulate( const boost::posix_time::ptime& time, const vector_type& mean ) = 0;
-
-		virtual background_model* clone() const = 0;
-};
-
-class background_daytime_model: public background_model {
-	private:
-		size_type   m_size;
-	public:
-		background_daytime_model( const boost::posix_time::ptime& time );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_twilight_model: public background_model {
-	public:
-		typedef boost::function1<value_type, const boost::posix_time::ptime&> sun_altitude_function_type;
-	private:
-		size_type   m_size;
-		vector_type m_background;
-		sun_altitude_function_type m_sun_altitude;
-	private:
-		value_type extrapolate( const boost::posix_time::ptime& timef ) const;
-	public:
-		background_twilight_model( const boost::posix_time::ptime& time, const sun_altitude_function_type& sun_altitude );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_night_model: public background_model {
-	private:
-		size_type   m_size;
-		vector_type m_background;
-	public:
-		background_night_model( const boost::posix_time::ptime& time );
-
-		size_type size() const;
-		vector_type background( const boost::posix_time::ptime& time ) const;
-		void accumulate( const boost::posix_time::ptime& time, const vector_type& mean );
-
-		background_model* clone() const;
-};
-
-class background_nonpoisson_filter {
-	public:
-		typedef background_factory::value_type  value_type;
-		typedef background_factory::vector_type vector_type;
-		typedef background_factory::size_type   size_type;
-
-		typedef background_factory::flux_type   flux_type;
-		typedef background_factory::vars_type   vars_type;
-	private:
-		value_type m_threshold;
-		value_type m_nonpoissonD;
-		value_type m_numexp;
-	public:
-		background_nonpoisson_filter( const value_type& threshold, const value_type& nonpoissonD, const value_type& numexp );
-		~background_nonpoisson_filter();
+typedef boost::posix_time::ptime   time_type;
+typedef double value_type;
+typedef vector<value_type>         vector_type;
+typedef vector_type::size_type     size_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
+
+typedef vector_type           flux_type;
+typedef symmetric_matrix_type vars_type;
+
+namespace model {
+	class abstract;
+}
+
+class factory {
+public:
+	typedef model::abstract              model_type;
+	typedef atmos::background::time_type time_type;
+	typedef atmos::background::flux_type flux_type;
+	
+	typedef boost::function1<model_type*, const time_type&> select_model_function_type;
+private:
+	std::auto_ptr< model_type > m_curmodel;
+	select_model_function_type  m_select_model;
+public:
+	factory(const select_model_function_type& select_model);
+	factory(const factory& factory);
+	void swap(factory& factory) throw();
+	factory& operator=(const factory& factory);
+	~factory();
+
+	flux_type background( const time_type& time ) const;
+
+	model_type* create_new_model( const time_type& time ) const;
+	model_type* release_current_model();
+	void handle_new_model( model_type* model );
+};
+
+namespace model {
+	
+class abstract {
+public:
+	typedef atmos::background::size_type size_type;
+	typedef atmos::background::time_type time_type;
+	typedef atmos::background::flux_type flux_type;
+private:
+	time_type m_inittime;
+public:
+	abstract( const time_type& time );
+	virtual ~abstract() = 0;
+
+	const time_type& inittime() const { return m_inittime; }

atmos-2.96.9.tar.gz/src/config.h.in -> atmos-2.97.1.tar.gz/src/config.h.in Changed

@@ -29,8 +29,9 @@
 /* Define to 1 if you have the <boost/io/ios_state.hpp> header file. */
 #undef HAVE_BOOST_IO_IOS_STATE_HPP
 
-/* Define to 1 if you have the <boost/iterator.hpp> header file. */
-#undef HAVE_BOOST_ITERATOR_HPP
+/* Define to 1 if you have the <boost/iterator/transform_iterator.hpp> header
+   file. */
+#undef HAVE_BOOST_ITERATOR_TRANSFORM_ITERATOR_HPP
 
 /* Define to 1 if you have the <boost/lexical_cast.hpp> header file. */
 #undef HAVE_BOOST_LEXICAL_CAST_HPP
@@ -38,8 +39,19 @@
 /* Define to 1 if you have the <boost/math/tools/roots.hpp> header file. */
 #undef HAVE_BOOST_MATH_TOOLS_ROOTS_HPP
 
-/* Define to 1 if you have the <boost/numeric/ublas/io.hpp> header file. */
-#undef HAVE_BOOST_NUMERIC_UBLAS_IO_HPP
+/* Define to 1 if you have the <boost/mem_fn.hpp> header file. */
+#undef HAVE_BOOST_MEM_FN_HPP
+
+/* Define to 1 if you have the <boost/mpl/if.hpp> header file. */
+#undef HAVE_BOOST_MPL_IF_HPP
+
+/* Define to 1 if you have the <boost/numeric/ublas/banded.hpp> header file.
+   */
+#undef HAVE_BOOST_NUMERIC_UBLAS_BANDED_HPP
+
+/* Define to 1 if you have the <boost/numeric/ublas/matrix_expression.hpp>
+   header file. */
+#undef HAVE_BOOST_NUMERIC_UBLAS_MATRIX_EXPRESSION_HPP
 
 /* Define to 1 if you have the <boost/numeric/ublas/matrix.hpp> header file.
    */
@@ -57,6 +69,10 @@
    file. */
 #undef HAVE_BOOST_NUMERIC_UBLAS_TRIANGULAR_HPP
 
+/* Define to 1 if you have the <boost/numeric/ublas/vector_expression.hpp>
+   header file. */
+#undef HAVE_BOOST_NUMERIC_UBLAS_VECTOR_EXPRESSION_HPP
+
 /* Define to 1 if you have the <boost/numeric/ublas/vector.hpp> header file.
    */
 #undef HAVE_BOOST_NUMERIC_UBLAS_VECTOR_HPP
@@ -71,6 +87,10 @@
 /* Define to 1 if you have the <boost/tuple/tuple.hpp> header file. */
 #undef HAVE_BOOST_TUPLE_TUPLE_HPP
 
+/* Define to 1 if you have the <boost/type_traits/is_complex.hpp> header file.
+   */
+#undef HAVE_BOOST_TYPE_TRAITS_IS_COMPLEX_HPP
+
 /* Define to 1 if you have the <gsl/gsl_errno.h> header file. */
 #undef HAVE_GSL_GSL_ERRNO_H
 
@@ -118,9 +138,6 @@
 /* Define to 1 if you have the <sys/types.h> header file. */
 #undef HAVE_SYS_TYPES_H
 
-/* Define to 1 if you have the <tr1/tuple> header file. */
-#undef HAVE_TR1_TUPLE
-
 /* Define to 1 if you have the <unistd.h> header file. */
 #undef HAVE_UNISTD_H

atmos-2.97.1.tar.gz/src/covariance_decomposition.h Added

@@ -0,0 +1,137 @@
+/*
+    $Id: covariance_decomposition.h 182 2011-01-31 11:26:03Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _COVARIANCE_DECOMPOSITION_H
+#define _COVARIANCE_DECOMPOSITION_H
+
+#include <cmath>
+
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/tuple/tuple.hpp>
+
+#include <lsp/cholesky_decomposition.h>
+
+#include "apply_to_all.h"
+
+namespace utils {
+namespace ublas {
+
+using namespace boost::numeric::ublas;
+
+namespace detail {
+	template<class M> boost::tuple< banded_matrix< typename M::value_type >, symmetric_matrix< typename M::value_type > > correlation_matrix( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		symmetric_matrix< value_type > cor(cov);
+		banded_matrix< value_type > sigma = utils::ublas::sqrt( banded_adaptor< symmetric_matrix< value_type > >(cor) );
+		for( typename banded_matrix< value_type >::size_type i = 0; i < sigma.size1(); ++i ) {
+			row(cor,i)    /= sigma(i,i);
+			column(cor,i) /= sigma(i,i);
+		}
+
+		return boost::make_tuple( sigma, cor );
+	}
+	template<class M> triangular_matrix< typename M::value_type > sqrt_matrix_decomposition( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		symmetric_matrix< value_type > decomp(cov);
+		lsp::cholesky_decomposition( decomp );
+		triangular_adaptor< symmetric_matrix< value_type > > result(decomp);
+		
+		return result;
+	}
+	template<class M> triangular_matrix< typename M::value_type > inv_sqrt_matrix_decomposition( const M& cov ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		triangular_matrix< value_type > triang( sqrt_matrix_decomposition(cov) );
+		matrix< value_type > inv = solve( triang, identity_matrix< value_type >( triang.size1() ), lower_tag() );
+
+		return triangular_adaptor< matrix< value_type > >(inv);
+	}
+	template<class M> triangular_matrix< typename M::value_type > aetkin_weight_matrix( const M& cov, const typename M::value_type& cond_threshold = 10000 ) {
+		typedef M matrix_type;
+		typedef typename M::value_type value_type;
+		typedef typename M::size_type  size_type;
+
+		if( BOOST_UBLAS_SAME( cov.size1(), cov.size2() ) < 1 )
+			return triangular_matrix< value_type >(0,0);
+
+		for( size_type i = 0; i < BOOST_UBLAS_SAME( cov.size1(), cov.size2() ); ++i ) {
+			if( cov(i,i) <= value_type(0) )
+				throw std::logic_error("Non-positive diagonal value");
+		}
+
+		boost::tuple< banded_matrix< value_type >, symmetric_matrix< value_type > > sigma_cor = correlation_matrix(cov);
+		triangular_matrix< value_type > weight = sqrt_matrix_decomposition( sigma_cor.template get<1>() );
+
+		value_type min = std::abs( weight(0,0) );
+		value_type max = std::abs( weight(0,0) );
+		for( typename triangular_matrix< value_type >::size_type i = 1; i < weight.size1(); ++i ) {
+			if( std::abs( weight(i,i) ) > max ) max = std::abs( weight(i,i) );
+			if( std::abs( weight(i,i) ) < min ) min = std::abs( weight(i,i) );
+		}
+
+		value_type cond = max / min;
+		if( !std::isfinite( cond ) || cond > cond_threshold ) {
+			throw std::logic_error( (boost::format("The condition number of covariance matrix is %1%") % cond).str() );
+		}
+
+		matrix< value_type > inv = solve( weight, identity_matrix< value_type >( weight.size1() ), lower_tag() );
+		inv = prod( inv, utils::ublas::inv( sigma_cor.template get<0>() ) );
+		return triangular_adaptor< matrix< value_type > >(inv);
+	}
+}
+
+template<class T> boost::tuple< banded_matrix<T>, symmetric_matrix<T> > correlation_matrix( const symmetric_matrix<T>& cov ) {
+	return detail::correlation_matrix( cov );
+}
+template<class T> triangular_matrix<T> sqrt_matrix_decomposition( const symmetric_matrix<T>& cov ) {
+	return detail::sqrt_matrix_decomposition( cov );
+}
+template<class T> triangular_matrix<T> inv_sqrt_matrix_decomposition( const symmetric_matrix<T>& cov ) {
+	return detail::inv_sqrt_matrix_decomposition( cov );
+}
+template<class T> triangular_matrix<T> aetkin_weight_matrix( const symmetric_matrix<T>& cov ) {
+	return detail::aetkin_weight_matrix( cov );
+}
+
+template<class M> boost::tuple< banded_matrix< typename M::value_type >, symmetric_matrix< typename M::value_type > > correlation_matrix( const symmetric_adaptor<M>& cov ) {
+	return detail::correlation_matrix( cov );
+}
+template<class M> triangular_matrix< typename M::value_type > sqrt_matrix_decomposition( const symmetric_adaptor<M>& cov ) {
+	return detail::sqrt_matrix_decomposition( cov );
+}
+template<class M> triangular_matrix< typename M::value_type > inv_sqrt_matrix_decomposition( const symmetric_adaptor<M>& cov ) {
+	return detail::inv_sqrt_matrix_decomposition( cov );
+}
+template<class M> triangular_matrix< typename M::value_type >  aetkin_weight_matrix( const symmetric_adaptor<M>& cov ) {
+	return detail::aetkin_weight_matrix( cov );
+}
+
+}}
+
+#endif // _COVARIANCE_DECOMPOSITION_H

atmos-2.96.9.tar.gz/src/dimm.cpp -> atmos-2.97.1.tar.gz/src/dimm.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: dimm.cpp 115 2010-03-14 08:48:08Z matwey $
+    $Id: dimm.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,13 +17,15 @@
 */
 
 #include <cmath>
-
-#include <boost/bind.hpp>
+#include <boost/mem_fn.hpp>
 
 #include "dimm.h"
 
-dimm::dimm_vars_type dimm::correction( const dimm_vars_type& vars ) {
-	dimm_vars_type ret( vars );
+namespace atmos {
+namespace dimm {
+
+vars_type correction( const vars_type& vars ) {
+	vars_type ret( vars );
 
 	ret(1) = ret(1)*ret(1) - ret(3)*ret(3); // CCD noise removal
 	value_type rho = ret(2)/ret(1);
@@ -33,48 +35,24 @@
 	return ret;
 }
 
-dimm::value_type dimm::turbulence( value_type scale, value_type Knumber, value_type variance ) {
+value_type turbulence( value_type scale, value_type Knumber, value_type variance ) {
 	value_type ret = variance;
 	ret *= scale * scale / 16.7 / Knumber;
 	return ret;
 }
 
-dimm::value_type dimm::seeing( value_type turbulence_power ) {
+value_type seeing( value_type turbulence_power ) {
 	value_type ret = 2.0e7 * std::pow( turbulence_power, 0.6 );
 	return ret;
 }
 
-const dimm::dimm_flux_type& dimm::left_flux( const dimm_sample& x ) {
-	return x.left_flux;
-}
-const dimm::dimm_flux_type& dimm::right_flux( const dimm_sample& x ) {
-	return x.right_flux;
-}
-const dimm::dimm_vars_type& dimm::xvars( const dimm_sample& x ) {
-	return x.xvars;
-}
-const dimm::dimm_vars_type& dimm::yvars( const dimm_sample& x ) {
-	return x.yvars;
-}
-
-dimm_storage::dimm_storage( size_type reserve ):
-	table_storage< dimm::dimm_sample >( reserve ),
-	m_left_flux_column( *this, boost::bind(&dimm::left_flux,_1) ),
-	m_right_flux_column( *this, boost::bind(&dimm::right_flux,_1) ),
-	m_xvars_column( *this, boost::bind(&dimm::xvars,_1) ),
-	m_yvars_column( *this, boost::bind(&dimm::yvars,_1) )
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_left_flux_column( *this, boost::mem_fn(&detail::sample::left_flux) ),
+	m_right_flux_column( *this, boost::mem_fn(&detail::sample::right_flux) ),
+	m_xvars_column( *this, boost::mem_fn(&detail::sample::xvars) ),
+	m_yvars_column( *this, boost::mem_fn(&detail::sample::yvars) )
 	{
 }
 
-void dimm_storage::correct_low_vars( ) {
-	double meanx = 0;
-	double meany = 0;
-	for( const_iterator it = begin(); it != end(); ++it ) {
-		meanx += it->xvars(0) / size();
-		meany += it->yvars(0) / size();
-	}
-	for( iterator it = begin(); it != end(); ++it ) {
-		it->xvars(1) += ( it->xvars(0) - meanx ) * ( it->xvars(0) - meanx );
-		it->yvars(1) += ( it->yvars(0) - meany ) * ( it->yvars(0) - meany );
-	}
-}
+}}

atmos-2.96.9.tar.gz/src/dimm.h -> atmos-2.97.1.tar.gz/src/dimm.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: dimm.h 115 2010-03-14 08:48:08Z matwey $
+    $Id: dimm.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,71 +19,99 @@
 #ifndef _DIMM_H
 #define _DIMM_H
 
-#include "table_storage.h"
-
 #include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
+#include "table_storage.h"
+
+namespace atmos {
 namespace dimm {
-	using namespace boost::numeric::ublas;
 
-	typedef boost::posix_time::ptime time_type;
-	
-	typedef double                 value_type;
-	typedef vector<value_type>     vector_type;
-	typedef vector_type::size_type size_type;
-
-	typedef vector_type dimm_flux_type;
-	typedef vector_type dimm_vars_type;
-
-	dimm_vars_type correction( const dimm_vars_type& vars );
-	value_type turbulence( value_type scale, value_type Knumber, value_type variance );
-	value_type seeing( value_type turbulence_power );
-
-	struct dimm_sample {
-		time_type      timestamp;
-		unsigned int   framenum;
-		dimm_flux_type left_flux;
-		dimm_flux_type right_flux;
-		dimm_vars_type xvars;
-		dimm_vars_type yvars;
-
-		dimm_sample( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ):
-			timestamp(time), framenum(framenum), left_flux(left), right_flux(right), xvars(xvars), yvars(yvars) {}
-	};
-
-	const dimm_flux_type& left_flux( const dimm_sample& x );
-	const dimm_flux_type& right_flux( const dimm_sample& x );
-	const dimm_vars_type& xvars( const dimm_sample& x );
-	const dimm_vars_type& yvars( const dimm_sample& x );
+using namespace boost::numeric::ublas;
+
+typedef boost::posix_time::ptime time_type;
+typedef double                   value_type;
+typedef vector<value_type>       vector_type;
+typedef vector_type::size_type   size_type;
+
+typedef vector_type  flux_type;
+typedef vector_type  vars_type;
+typedef unsigned int franenum_type;
+
+vars_type  correction( const vars_type& vars );
+value_type turbulence( value_type scale, value_type Knumber, value_type variance );
+value_type seeing( value_type turbulence_power );
+
+template<class Iterator, class T> void correct_low_vars( Iterator begin, Iterator end, const T& mean ) {
+	typedef Iterator iterator_type;
+	typedef T value_type;
+
+	for( iterator_type it = begin; it != end; ++it ) {
+		(*it)(1) += ((*it)(0)-mean)*((*it)(0)-mean);
+	}
+}
+
+namespace detail {
+
+struct sample {
+	typedef atmos::dimm::time_type     timestamp_type;
+	typedef atmos::dimm::franenum_type framenum_type;
+	typedef atmos::dimm::flux_type     flux_type;
+	typedef atmos::dimm::vars_type     vars_type;
+		
+	timestamp_type timestamp;
+	framenum_type  framenum;
+	flux_type      left_flux;
+	flux_type      right_flux;
+	vars_type      xvars;
+	vars_type      yvars;
+
+	sample( const timestamp_type& time,
+		     const framenum_type& framenum,
+		     const flux_type& left,
+		     const flux_type& right,
+		     const vars_type& xvars,
+		     const vars_type& yvars ):
+		timestamp(time), framenum(framenum), left_flux(left), right_flux(right), xvars(xvars), yvars(yvars) {}
 };
+}
 
-class dimm_storage:
-	public table_storage< dimm::dimm_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_flux_type, dimm::dimm_sample > > left_flux_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_flux_type, dimm::dimm_sample > > right_flux_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_vars_type, dimm::dimm_sample > > xvars_column_type;
-	typedef table_storage_column< dimm::dimm_sample, boost::function1< dimm::dimm_vars_type, dimm::dimm_sample > > yvars_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::flux_type > left_flux_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::flux_type > right_flux_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::vars_type > xvars_column_type;
+	typedef table_storage_column_proxy< self_type, detail::sample::vars_type > yvars_column_type;
 private:
 	left_flux_column_type  m_left_flux_column;
 	right_flux_column_type m_right_flux_column;
 	xvars_column_type      m_xvars_column;
 	yvars_column_type      m_yvars_column;
 public:
-	dimm_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const left_flux_column_type&  left_flux()  const { return m_left_flux_column; }
 	const right_flux_column_type& right_flux() const { return m_right_flux_column; }
 	const xvars_column_type&      xvars()      const { return m_xvars_column; }
 	const yvars_column_type&      yvars()      const { return m_yvars_column; }
+	left_flux_column_type&  left_flux()  { return m_left_flux_column; }
+	right_flux_column_type& right_flux() { return m_right_flux_column; }
+	xvars_column_type&      xvars()      { return m_xvars_column; }
+	yvars_column_type&      yvars()      { return m_yvars_column; }
+
+	void push_back( const time_type& time,
+			const franenum_type& framenum,
+			const flux_type& left,
+			const flux_type& right,
+			const vars_type& xvars,
+			const vars_type& yvars ) {
 
-	void push_back( const dimm::time_type& time, unsigned int framenum, const dimm::dimm_flux_type& left, const dimm::dimm_flux_type& right, const dimm::dimm_vars_type& xvars, const dimm::dimm_flux_type& yvars ) {
-		// FIXME: dimm_sample_ref
-		table_storage< dimm::dimm_sample >::push_back( dimm::dimm_sample( time, framenum, left, right, xvars, yvars ) );
+		table_storage< detail::sample >::push_back( detail::sample( time, framenum, left, right, xvars, yvars ) );
 	}
-
-	void correct_low_vars();
 };
 
+}}
+
 #endif // _DIMM_H

atmos-2.96.9.tar.gz/src/fft.h -> atmos-2.97.1.tar.gz/src/fft.h Changed

atmos-2.96.9.tar.gz/src/file_input.cpp -> atmos-2.97.1.tar.gz/src/file_input.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_input.cpp 142 2010-04-03 13:50:39Z matwey $
+    $Id: file_input.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,34 +16,34 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <cassert>
-#include <stdexcept>
-#include <sstream>
 #include <iostream>
+#include <sstream>
+#include <stdexcept>
 
-#include <ctime>
-
-#include <boost/io/ios_state.hpp>
-#include <boost/numeric/ublas/vector.hpp>
-
+#include <boost/algorithm/string.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
 
 #include "atmos.h"
 #include "file_input.h"
+#include "io.h"
+
+namespace atmos {
+namespace input {
 
-file_input::file_input( std::istream& stm ):
+file::file( std::istream& stm ):
 	m_linenum(0),
 	m_stm(stm),
 	m_cached( false ),
 	m_locale( m_stm.getloc(), new boost::posix_time::time_input_facet("%Y-%m-%d %H:%M:%S") ) {
 }
-file_input::~file_input() {
+file::~file() {
 }
-bool file_input::unshift() throw() {
+
+bool file::unshift() throw() {
 	m_cached = true;
 	return true;
 }
-bool file_input::shift() {
+bool file::shift() {
 	struct linenum_guard {
 			linenum_guard( size_type& linenum ): m_linenum(linenum) {}
 			~linenum_guard() { m_linenum++; }
@@ -65,16 +65,15 @@
 	
 	return false;
 }
-file_input::size_type file_input::line() const throw() {
+size_type file::line() const throw() {
 	return m_linenum;
 }
-//FIXME: use boost string_algo
-std::string file_input::strip_spaces( const std::string& str ) {
+std::string file::strip_spaces( const std::string& str ) {
 	std::string ret(str);
 	boost::trim_if( ret, boost::is_space() );
 	return ret;
 }
-void file_input::parse_line( size_type linenum, const std::string& line ) {
+void file::parse_line( size_type linenum, const std::string& line ) {
 	if( line.size() == 0 )
 		return; // skip empty line
 
@@ -96,10 +95,10 @@
 		throw failure("Bad line format");
 	}
 }
-void file_input::parse_comment_line( std::istringstream& stm ) {
+void file::parse_comment_line( std::istringstream& stm ) {
 	on_comment( stm.str() );
 }
-void file_input::parse_param_line( std::istringstream& stm ) {
+void file::parse_param_line( std::istringstream& stm ) {
 	static time_type timestamp; // NOTE: thread unsafe
 	std::string key,value;
 	
@@ -109,7 +108,7 @@
 
 	on_parameter(timestamp, strip_spaces(key), strip_spaces(value));
 }
-void file_input::parse_object_line( std::istringstream& stm ) {
+void file::parse_object_line( std::istringstream& stm ) {
 	static time_type   timestamp; // NOTE: thread unsafe
 	static object_type obj;       // NOTE: thread unsafe
 
@@ -117,7 +116,7 @@
 
 	on_object( timestamp, obj );
 }
-void file_input::parse_mode_line( std::istringstream& stm ) {
+void file::parse_mode_line( std::istringstream& stm ) {
 	static time_type timestamp; // NOTE: thread unsafe
 	
 	stm >> timestamp;
@@ -132,7 +131,7 @@
 		throw failure("Unknown mode");
 	}
 }
-void file_input::parse_measurement_line( std::istringstream& stm ) {
+void file::parse_measurement_line( std::istringstream& stm ) {
 	static time_type timestamp;               // NOTE: thread unsafe
 	static flux_type flux( atmos::channels ); // NOTE: thread unsafe
 	static vars_type vars( atmos::channels, atmos::channels ); // NOTE: thread unsafe
@@ -147,7 +146,7 @@
 
 	on_measurements( timestamp, flux, vars, rho1, rho2 );
 }
-void file_input::parse_dimm_line( std::istringstream& stm ) {
+void file::parse_dimm_line( std::istringstream& stm ) {
 	static time_type timestamp;  // NOTE: thread unsafe
 	static unsigned int framenum;
 	static dimm_flux_type left_flux( 3 );  // NOTE: thread unsafe
@@ -169,7 +168,10 @@
 	on_dimm( timestamp, framenum, left_flux, right_flux, xvars, yvars );
 }
 
-file_input::failure::failure( const std::string& what):
+file::failure::failure( const std::string& what):
 	std::runtime_error( what ) {
 }
 
+}}
+
+

atmos-2.96.9.tar.gz/src/file_input.h -> atmos-2.97.1.tar.gz/src/file_input.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_input.h 142 2010-04-03 13:50:39Z matwey $
+    $Id: file_input.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -23,37 +23,49 @@
 #include <stdexcept>
 #include <string>
 
-#include <boost/algorithm/string.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
-
-#include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
-#include "io.h"
+#include "atmos.h" // size_type
+#include "typesdef.h" // object
 
-using namespace boost::numeric;
+namespace atmos {
+namespace input {
 
-class file_input {
-public:
-	typedef size_t size_type;
-	
-	typedef boost::posix_time::ptime time_type;
-	typedef ublas::vector<double>    flux_type;
-	typedef ublas::symmetric_matrix<double> vars_type;
-	typedef ublas::symmetric_matrix<double> rho1_type;
-	typedef ublas::vector<double>    rho2_type;
-	typedef object                   object_type;
+using namespace boost::numeric::ublas;
+
+typedef boost::posix_time::ptime time_type;
+typedef double               value_type;
+typedef vector< value_type > vector_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
+
+typedef vector_type flux_type;
+typedef symmetric_matrix_type vars_type;
+typedef symmetric_matrix_type rho1_type;
+typedef vector_type rho2_type;
+typedef vector_type dimm_flux_type;
+typedef vector_type dimm_vars_type;
+typedef object object_type;
 
-	typedef ublas::vector<double>    dimm_flux_type;
-	typedef ublas::vector<double>    dimm_vars_type;
+class file {
+public:
+	typedef atmos::size_type        size_type;
+	typedef atmos::input::time_type time_type;
+	typedef atmos::input::flux_type flux_type;
+	typedef atmos::input::vars_type vars_type;
+	typedef atmos::input::rho1_type rho1_type;
+	typedef atmos::input::rho2_type rho2_type;
+	typedef atmos::input::dimm_flux_type dimm_flux_type;
+	typedef atmos::input::dimm_vars_type dimm_vars_type;
+	typedef atmos::input::object_type object_type;
 private:
 	size_type m_linenum;
 	std::istream& m_stm;
 	std::string m_line;
 	bool m_cached;
 	std::locale m_locale;
-	private:
-	
+private:
 	static std::string strip_spaces( const std::string& str );
 
 	void parse_line( size_type linenum, const std::string& line );
@@ -63,7 +75,6 @@
 	void parse_mode_line( std::istringstream& stm );
 	void parse_measurement_line( std::istringstream& stm );
 	void parse_dimm_line( std::istringstream& stm );
-
 protected:
 	bool unshift() throw();
 public:
@@ -73,11 +84,10 @@
 		failure( const std::string& what );
 	};
 public:
-	file_input( std::istream& stm );
-	virtual ~file_input();
+	file( std::istream& stm );
+	virtual ~file();
 
 	bool shift();
-
 	size_type line() const throw() ;
 public:
 	virtual void on_comment( const std::string& line ) = 0;
@@ -90,4 +100,6 @@
 	virtual void on_dimm( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) = 0;
 };
 
+}}
+
 #endif // _FILE_INPUT_H

atmos-2.96.9.tar.gz/src/file_output.cpp -> atmos-2.97.1.tar.gz/src/file_output.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.cpp 108 2010-03-07 11:40:13Z matwey $
+    $Id: file_output.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,30 +16,36 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#include <boost/date_time/posix_time/posix_time.hpp>
+#include <boost/format.hpp>
 #include <boost/io/ios_state.hpp>
 
+#include "atmos.h"
 #include "file_output.h"
-#include "integrals.h"
+#include "io.h"
 
-file_output::file_output( std::ostream& stm ): m_stm(stm) {
+namespace atmos {
+namespace output {
+
+file::file( std::ostream& stm ): m_stm(stm) {
 	boost::posix_time::time_facet* facet(new boost::posix_time::time_facet("%Y-%m-%d %H:%M:%S"));
 	stm.imbue(std::locale(stm.getloc(), facet));
 }
-file_output::~file_output() {
+file::~file() {
 
 }
-void file_output::comment( const std::string& comment ) {
+void file::comment( const std::string& comment ) {
 	m_stm << "# " << comment << std::endl;
 }
-void file_output::normal( const time_type& time, double mz ) {
+void file::normal( const time_type& time, double mz ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "M " << time << " Normal " << std::fixed << std::setprecision(3) << mz << std::endl;
 }
-void file_output::background( const time_type& time ) {
+void file::background( const time_type& time ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "M " << time << " Background measurements" << std::endl;
 }
-void file_output::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 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_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 time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err ) {
 	static datum skipped;
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "A " << time << " "
@@ -50,7 +56,7 @@
 		<< std::scientific << datum( moment_2, moment_2_err/moment_2 ) << " "
 		<< std::fixed << datum( time_constant, time_constant_err/time_constant ) << std::endl;
 }
-void file_output::profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const ublas::vector<double>& cn2, const ublas::vector<double>& ecn2, const ublas::vector<double>& grid ){
+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() );
 	assert( cn2.size() == grid.size() );
 
@@ -58,27 +64,27 @@
 	m_stm << "T " << time << " " << method << " "
 		<< std::setw(2) << cn2.size() << " " << std::setprecision(2)
 		<< std::fixed << chi2 << " " << cn2seeing;
-	for( ublas::vector<double>::size_type i = 0; i < cn2.size(); i++)
+	for( vector_type::size_type i = 0; i < cn2.size(); i++)
 		m_stm << " " << std::fixed << grid(i)
 			<< std::scientific  << " " << cn2(i) << " " << ecn2(i);
 	m_stm << std::endl;
 }
-void file_output::residual( const time_type& time, const std::string& method, const ublas::vector<double>& o_csi ) {
+void file::residual( const time_type& time, const std::string& method, const vector_type& o_csi ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "R " << time << " " << method << " "
 		<< std::setprecision(3) << std::fixed;
-	for( ublas::vector<double>::size_type i = 0; i < o_csi.size(); i++ )
+	for( vector_type::size_type i = 0; i < o_csi.size(); i++ )
 		m_stm << "  " << o_csi(i);
 	m_stm << std::endl;
 }
-void file_output::fluxes( const time_type& time, const ublas::vector<double>& result_flux ) {
+void file::fluxes( const time_type& time, const vector_type& result_flux ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	assert( result_flux.size() >= 4 );
 	m_stm << "F " << time << " " << std::fixed << datum( result_flux(0) ) << ' ' << datum( result_flux(1) ) << ' ' << datum( result_flux(2) ) << ' ' << datum( result_flux(3) ) << std::endl;
 }
-void file_output::fluxes( const time_type& time, const ublas::vector<double>& result_flux, const ublas::vector<double>& result_flux_err ) {
+void file::fluxes( const time_type& time, const vector_type& result_flux, const vector_type& result_flux_err ) {
 	boost::io::ios_all_saver asaver( m_stm );
-	assert( result_flux.size() >= 4 );
+	assert( result_flux.size() == atmos::channels );
 	assert( result_flux.size() == result_flux_err.size() );
 	m_stm << "F " << time << " " << std::fixed
 	<< datum( result_flux(0), result_flux_err(0)/result_flux(0) ) << ' '
@@ -86,12 +92,12 @@
 	<< datum( result_flux(2), result_flux_err(2)/result_flux(2) ) << ' '
 	<< datum( result_flux(3), result_flux_err(3)/result_flux(3) ) << std::endl;
 }
-void file_output::indices( const time_type& time, const ublas::symmetric_matrix<double>& result_indx, const ublas::symmetric_matrix<double>& result_indx_err, const ublas::vector<double>& result_desi, const ublas::vector<double>& result_desi_err ) {
+void file::indices( const time_type& time, const symmetric_matrix_type& result_indx, const symmetric_matrix_type& result_indx_err, const vector_type& result_desi, const vector_type& result_desi_err ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	ios_indic_saver isaver( m_stm );
-	assert( result_desi.size() >= 4 );
-	assert( result_indx.size1() >= 4 );
-	assert( result_indx.size2() >= 4 );
+	assert( result_desi.size() == atmos::channels );
+	assert( result_indx.size1() == atmos::channels );
+	assert( result_indx.size2() == atmos::channels );
 	assert( result_desi.size() == result_desi_err.size() );
 	assert( result_indx.size1() == result_indx_err.size1() );
 	assert( result_indx.size1() == result_indx_err.size2() );
@@ -102,11 +108,13 @@
 		datum( result_indx(2,3), result_indx_err(2,3)/result_indx(2,3) ) <<
 		datum( result_desi(0), result_desi_err(0)/result_desi(0) ) << datum( result_desi(1), result_desi_err(1)/result_desi(1) ) << datum( result_desi(2), result_desi_err(2)/result_desi(2) ) << datum( result_desi(3), result_desi_err(3)/result_desi(3) ) << std::endl;
 }
-void file_output::parameter( const time_type& time, const std::string& key, const std::string& value ) {
+void file::parameter( const time_type& time, const std::string& key, const std::string& value ) {
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "P " << time << " " << key << " = " << value << std::endl;
 }
-void file_output::objectinfo( const time_type& time, const object_type& obj ){
+void file::objectinfo( const time_type& time, const object_type& obj ){
 	boost::io::ios_all_saver asaver( m_stm );
 	m_stm << "O " << time << " " << obj << std::endl;
 };
+
+}}

atmos-2.96.9.tar.gz/src/file_output.h -> atmos-2.97.1.tar.gz/src/file_output.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: file_output.h 140 2010-04-02 17:13:12Z matwey $
+    $Id: file_output.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,42 +22,46 @@
 #include <iosfwd>
 #include <string>
 
-#include <boost/format.hpp>
-#include <boost/io/ios_state.hpp>
+#include <boost/date_time/posix_time/posix_time.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 
-#include <boost/date_time/posix_time/posix_time.hpp>
+#include "typesdef.h" // object
 
-#include "io.h"
+namespace atmos {
+namespace output {
 
-using namespace boost::numeric;
+using namespace boost::numeric::ublas;
 
-class file_output {
-public:
-	typedef size_t size_type;
+typedef boost::posix_time::ptime time_type;
+typedef double               value_type;
+typedef vector< value_type > vector_type;
+typedef symmetric_matrix< value_type > symmetric_matrix_type;
 
-	typedef boost::posix_time::ptime time_type;
-	typedef object                   object_type;
+typedef object object_type;
+
+class file {
 private:
 	std::ostream& m_stm;
 public:
-	file_output( std::ostream& stm );
-	~file_output();
+	file( std::ostream& stm );
+	~file();
 
 	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 time_constant, double time_constant_err, double moment_0, double moment_0_err, double moment_2, double moment_2_err );
-	void profile( const time_type& time, const std::string& method, double cn2seeing, double chi2, const ublas::vector<double>& cn2, const ublas::vector<double>& ecn2, const ublas::vector<double>& grid );
-	void residual( const time_type& time, const std::string& method, const ublas::vector<double>& o_csi );
-	void fluxes( const time_type& time, const ublas::vector<double>& result_flux );
-	void fluxes( const time_type& time, const ublas::vector<double>& result_flux, const ublas::vector<double>& result_flux_err );
-	void indices( const time_type& time, const ublas::symmetric_matrix<double>& result_indx, const ublas::symmetric_matrix<double>& result_indx_err, const ublas::vector<double>& result_desi, const ublas::vector<double>& result_desi_err );
+	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 );
+	void fluxes( const time_type& time, const vector_type& result_flux, const vector_type& result_flux_err );
+	void indices( const time_type& time, const symmetric_matrix_type& result_indx, const symmetric_matrix_type& result_indx_err, const vector_type& result_desi, const vector_type& result_desi_err );
 	void parameter( const time_type& time, const std::string& key, const std::string& value );
 	void objectinfo( const time_type& time, const object_type& obj );
 
 	std::ostream& operator()() { return m_stm; }
 };
 
+}}
+
 #endif // _FILE_OUTPUT_H

atmos-2.96.9.tar.gz/src/gsl_wrap.h -> atmos-2.97.1.tar.gz/src/gsl_wrap.h Changed

atmos-2.96.9.tar.gz/src/indices.cpp -> atmos-2.97.1.tar.gz/src/indices.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: indices.cpp 147 2010-04-17 09:25:25Z matwey $
+    $Id: indices.cpp 181 2011-01-31 10:07:58Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,28 +16,27 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <algorithm>
-#include <numeric>
+#include <stdexcept>
 
 #include <boost/format.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 
-#include "atmos.h"
 #include "indices.h"
 #include "io.h"
+#include "photometry.h"
 #include "utils.h"
 
-#include "photometry.h"
+namespace atmos {
+namespace indices {
 
-indices::flux_type indices::flux_nonlinearity_correction ( const flux_type& flux, const vector_type& deadtime ) {
+flux_type flux_nonlinearity_correction ( const flux_type& flux, const vector_type& deadtime ) {
 	return photometry::flux_nonlinearity_correction( flux, deadtime );
 }
-indices::flux_type indices::flux_background_correction( const flux_type& flux, const vector_type& scatter, const vector_type& background ) {
+flux_type flux_background_correction( const flux_type& flux, const vector_type& scatter, const vector_type& background ) {
 	return photometry::flux_scatter_correction( flux, scatter ) - background;
 }
-boost::tuple< indices::vars_type, indices::rho1_type, indices::rho2_type >
-	indices::vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson ) {
+boost::tuple< vars_type, rho1_type, rho2_type > vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson ) {
 	using namespace utils::ublas;
 
 	typedef matrix_vector_slice< vars_type > diagonal_type;
@@ -65,11 +64,10 @@
 	return ret;
 }
 
-boost::tuple< indices::indices_type, indices::indesis_type >
-	indices::calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+boost::tuple< indices_type, indesis_type > calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
 	using namespace utils::ublas;
 
-	typedef matrix_vector_slice< indices_type >    indices_diagonal_type;
+	typedef matrix_vector_slice< const indices_type > indices_diagonal_type;
 	typedef matrix_vector_slice< const rho1_type > rho1_diagonal_type;
 	typedef boost::tuple< indices_type, indesis_type> result_type;
 
@@ -82,29 +80,27 @@
 	
 	rho1_diagonal_type rho1_diag = diag( rho1 );
 
-	ret.get<0>() = element_div( vars*1.25 - rho1*0.25, outer_prod( flux, flux ) );
-
-	indices_diagonal_type variances = diag( ret.get<0>() );
-
-	ret.get<1>() = ( variances*3.0 + element_div( rho2 - rho1_diag * 4.0, element_prod(flux,flux) ) ) / 4.5;
-
-	ret.get<0>() = utils::ublas::log( xp1( ret.get<0>() ) );
-
-	triangle( ret.get<0>(), unit_lower() ) = triangle(
-			outer_prod( variances, scalar_vector< value_type >( variances.size(), 1.0 ) )
-			+ outer_prod( scalar_vector< value_type >( variances.size(), 1.0 ), variances )
-			- triangle( ret.get<0>(), unit_lower() ) * 2.0,
-		unit_lower() );
+	ret.get<0>() = utils::ublas::log( xp1( element_div( vars*1.17 - rho1*0.17, outer_prod( flux, flux ) ) ) );
+	ret.get<1>() = ( diag( element_div( vars, outer_prod( flux, flux ) ) )*3.0 + element_div( rho2 - rho1_diag * 4.0, element_prod(flux,flux) ) ) / 4.5;
 
 	return ret;
 }
 
-indices_flux_threshold_filter::indices_flux_threshold_filter( value_type threshold ):
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_flux_column( *this, boost::mem_fn(&detail::sample::flux) ),
+	m_indices_column( *this, boost::mem_fn(&detail::sample::indices) ),
+	m_indesis_column( *this, boost::mem_fn(&detail::sample::indesis) ) {
+}
+
+namespace filter {
+	
+flux_threshold::flux_threshold( value_type threshold ):
 	m_threshold( threshold ) {
 }
-indices_flux_threshold_filter::~indices_flux_threshold_filter() {
+flux_threshold::~flux_threshold() {
 }
-bool indices_flux_threshold_filter::operator() ( const flux_type& flux ) const {
+bool flux_threshold::operator() ( const flux_type& flux ) const {
 	if( *std::min_element( flux.begin(), flux.end() ) < m_threshold ) {
 		throw std::logic_error( (boost::format("Object flux %.2f is under the threshold %.2f") % flux % m_threshold).str() );
 		return true;
@@ -112,19 +108,5 @@
 	return false;
 }
 
-const indices::flux_type& indices::flux( const indices::indices_sample& x ) {
-	return x.flux;
-}
-const indices::indices_type& indices::indexes( const indices::indices_sample& x ) {
-	return x.indices;
-}
-const indices::indesis_type& indices::indesis( const indices::indices_sample& x ) {
-	return x.indesis;
-}
-
-indices_storage::indices_storage( size_type reserve ):
-	table_storage< indices::indices_sample >( reserve ),
-	m_flux_column( *this, boost::bind(&indices::flux,_1) ),
-	m_indices_column( *this, boost::bind(&indices::indexes,_1) ),
-	m_indesis_column( *this, boost::bind(&indices::indesis,_1) ) {
 }
+}}

atmos-2.96.9.tar.gz/src/indices.h -> atmos-2.97.1.tar.gz/src/indices.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: indices.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: indices.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,29 +19,23 @@
 #ifndef _INDICES_H
 #define _INDICES_H
 
-#include <stdexcept>
-#include <list>
-#include <vector>
-#include <iterator>
-
-#include <boost/bind.hpp>
 #include <boost/date_time/posix_time/posix_time.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/tuple/tuple.hpp>
 
 #include "table_storage.h"
 
+namespace atmos {
 namespace indices {
 	using namespace boost::numeric::ublas;
 
 	typedef boost::posix_time::ptime time_type;
-	
-	typedef double                value_type;
-	typedef vector< value_type >  vector_type;
-	typedef vector_type::size_type size_type;
+	typedef double                   value_type;
+	typedef vector< value_type >     vector_type;
+	typedef vector_type::size_type   size_type;
 	typedef symmetric_matrix< value_type > symmetric_matrix_type;
+	
 	typedef vector_type           flux_type;
 	typedef symmetric_matrix_type vars_type;
 	typedef symmetric_matrix_type rho1_type;
@@ -55,27 +49,33 @@
 		vars_correction( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2, const vector_type& deadtime, const vector_type& npoisson );
 	boost::tuple< indices_type, indesis_type > calculate( const flux_type& flux, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 );
 
-	struct indices_sample {
-		time_type    timestamp;
-		flux_type    flux;
-		indices_type indices;
-		indesis_type indesis;
-
-		indices_sample( const time_type& time, const flux_type& flux, const indices_type& indices, const indesis_type& indesis ):
+namespace detail {
+	struct sample {
+		typedef atmos::indices::time_type    timestamp_type;
+		typedef atmos::indices::flux_type    flux_type;
+		typedef atmos::indices::indices_type indices_type;
+		typedef atmos::indices::indesis_type indesis_type;
+
+		timestamp_type timestamp;
+		flux_type      flux;
+		indices_type   indices;
+		indesis_type   indesis;
+
+		sample( const timestamp_type& time,
+			const flux_type& flux,
+			const indices_type& indices,
+			const indesis_type& indesis ):
 			timestamp(time), flux(flux), indices(indices), indesis(indesis) {}
 	};
+}
 
-	const flux_type& flux( const indices_sample& x );
-	const indices_type& indexes( const indices_sample& x );
-	const indesis_type& indesis( const indices_sample& x );
-};
-
-class indices_storage:
-	public table_storage< indices::indices_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::flux_type, indices::indices_sample > >    flux_column_type;
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::indices_type, indices::indices_sample > > indices_column_type;
-	typedef table_storage_column< indices::indices_sample, boost::function1< indices::indesis_type, indices::indices_sample > > indesis_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::flux_type >    flux_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::indices_type > indices_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::indesis_type > indesis_column_type;
 
 private:
 	flux_column_type    m_flux_column;
@@ -83,28 +83,35 @@
 	indesis_column_type m_indesis_column;
 
 public:
-	indices_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const flux_column_type&    flux()    const { return m_flux_column; }
 	const indices_column_type& indices() const { return m_indices_column; }
 	const indesis_column_type& indesis() const { return m_indesis_column; }
 
-	void push_back( const indices::time_type& time, const indices::flux_type& flux, const indices::indices_type& indices, const indices::indesis_type& indesis ) {
-		table_storage< indices::indices_sample >::push_back( indices::indices_sample( time, flux, indices, indesis ) );
+	void push_back( const time_type& time,
+			const flux_type& flux,
+			const indices_type& indices,
+			const indesis_type& indesis ) {
+		
+		table_storage< detail::sample >::push_back( detail::sample( time, flux, indices, indesis ) );
 	}
 };
 
-class indices_flux_threshold_filter {
+namespace filter {
+	class flux_threshold {
 	public:
-		typedef indices::value_type value_type;
-		typedef indices::flux_type  flux_type;
+		typedef atmos::indices::flux_type  flux_type;
+		typedef atmos::indices::value_type value_type;
 	private:
 		value_type m_threshold;
 	public:
-		indices_flux_threshold_filter( value_type threshold );
-		~indices_flux_threshold_filter();
+		flux_threshold( value_type threshold );
+		~flux_threshold();
 
 		bool operator() ( const flux_type& flux ) const;
-};
+	};
+}
+}}
 
 #endif // _INDICES_H

atmos-2.96.9.tar.gz/src/integrals.cpp -> atmos-2.97.1.tar.gz/src/integrals.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: integrals.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,18 +17,24 @@
 */
 
 #include <cmath>
-#include <algorithm>
-#include <numeric>
-#include <boost/bind.hpp>
+
 #include <boost/format.hpp>
+#include <boost/mem_fn.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
+
 #include <lsp/singular_decomposition.h>
 
+#include "atmos.h"
 #include "apply_to_all.h"
-#include "indices.h"
 #include "integrals.h"
+#include "vector_generator.h"
+#include "utils.h"
+#include "weif.h"
+
+namespace atmos {
+namespace integrals {
 
 const double KM2M = 1000;
 const double POWSEE = 0 ;       // Power of Cn2 moment for seeing
@@ -64,11 +70,12 @@
 
 const double SPOW_TC  = -0.6;   // Power of DESI for atm. time constant
 
-integrals_factory::dcf_matrix_type integrals_factory::calculate( const std::string& spectral_type, const weight_type& awf ) {
+functor_factory::dcf_matrix_type functor_factory::calculate( const std::string& spectral_type, const weight_type& awf ) {
 	using namespace boost::numeric::ublas;
+	using namespace utils::algo;
 	using namespace utils::ublas;
 
-	const vector<double> grid( utils::ublas::vector_generator< lg_generator > ( lg_generator(0.0, 0.04, 0.25), 50 ) );
+	const vector<double> grid( 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;
@@ -81,7 +88,7 @@
 	assert( lower != grid.end() );
 	assert( higher != grid.rend() );
 
-	result_type res( npower, 12-2/*dimm*/ );
+	result_type res( npower, 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, awf) );
@@ -94,7 +101,7 @@
 		vector<double> alt    = utils::ublas::exp( log( h )*moment_power ) * scale;
 		vector<double> weight = utils::ublas::exp( log( h )*weight_power );
 
-		matrix<double> A( trans( subrange( wfv, 0, 10, lower.index(), higher.index()+1 ) ) );
+		matrix<double> A( trans( subrange( wfv, 0, atmos::indices_size, lower.index(), higher.index()+1 ) ) );
 		matrix<double> left( identity_matrix<double>( A.size1() ) );
 		matrix<double> right( identity_matrix<double>( A.size2() ) );
 
@@ -121,22 +128,22 @@
 
 	return res;
 }
-integrals_factory::integrals_factory() {
+functor_factory::functor_factory() {
 }
-integrals_factory::integrals_factory( const integrals_factory& factory ) {
+functor_factory::functor_factory( const functor_factory& factory ) {
 	m_storage = factory.m_storage;
 }
-void integrals_factory::swap( integrals_factory& factory ) throw() {
+void functor_factory::swap( functor_factory& factory ) throw() {
 	m_storage.swap( factory.m_storage );
 }
-integrals_factory& integrals_factory::operator=( const integrals_factory& factory ) {
-	integrals_factory tmp( factory );
+functor_factory& functor_factory::operator=( const functor_factory& factory ) {
+	functor_factory tmp( factory );
 	tmp.swap( *this );
 	return *this;
 }
-integrals_factory::~integrals_factory() {
+functor_factory::~functor_factory() {
 }
-integrals_factory::functor_type integrals_factory::get_functor( const std::string& spectral_type, const weight_type& awf ) {
+functor_factory::functor_type functor_factory::operator() ( const std::string& spectral_type, const weight_type& awf ) {
 	storage_type::const_iterator it = m_storage.find( spectral_type );
 	if( it != m_storage.end() )
 		return functor_type( it->second );
@@ -146,12 +153,15 @@
 	return functor_type( ret.first->second );
 }
 
-integrals_functor::result_type integrals_functor::operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const {
-	assert( indices.size1() == 4 );
-	assert( indices.size2() == 4 );
-	assert( indesis.size() == 4 );
+functor::result_type functor::operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const {
+	typedef vector< double >       vector_type;
+	typedef vector_type::size_type size_type;
+	
+	assert( indices.size1() == atmos::channels );
+	assert( indices.size2() == atmos::channels );
+	assert( indesis.size()  == atmos::channels );
 
-	vector< double > aints(10);
+	vector_type aints( atmos::indices_size );
 	aints(0) = indices(0,0); aints(1) = indices(1,1); aints(2) = indices(2,2); aints(3) = indices(3,3);
 	aints(4) = indices(0,1); aints(5) = indices(0,2); aints(6) = indices(0,3);
 	aints(7) = indices(1,2); aints(8) = indices(1,3); aints(9) = indices(2,3);
@@ -159,52 +169,31 @@
 	aints = prod( m_dcf, aints );
 
 	assert( npower == aints.size() );
-	for( vector< double >::size_type i = 0; i < npower; ++i ){
+	for( size_type i = 0; i < npower; ++i ){
 		aints(i) = aints(i) / std::pow( mz, 1.0 + powersi );
 	}
 
-	result_type ret;
-
-	ret.free_seeing       = KS * LEFF_SEE * std::pow( aints(0), JPOW_SEE );
-	ret.isoplanatic_angle = KP * LEFF_ISP * std::pow( aints(1), JPOW_ISP );
-	ret.altitude_eff      = 0.001 * std::pow( aints(1) / aints(2), 1.0/( POWISP-POWEFF ) );
-	ret.time_constant     = KT * std::pow( indesis(0), SPOW_TC );
-	ret.moment_0          = aints(0);
-	ret.moment_2          = aints(3);
-
-	return ret;
-}
-
-double integrals::free_seeing( const integrals_sample& x ) {
-	return x.free_seeing;
-}
-double integrals::isoplanatic_angle( const integrals_sample& x ) {
-	return x.isoplanatic_angle;
-}
-double integrals::altitude_eff( const integrals_sample& x ) {
-	return x.altitude_eff;
-}
-double integrals::time_constant( const integrals_sample& x ) {
-	return x.time_constant;
-}
-double integrals::moment_0( const integrals_sample& x ) {
-	return x.moment_0;
-}
-double integrals::moment_2( const integrals_sample& x ) {
-	return x.moment_2;
-}
-
-integrals_storage::integrals_storage( size_type reserve ):
-	table_storage< integrals::integrals_sample >( reserve ),
-	m_free_seeing_column( *this, boost::bind(&integrals::free_seeing,_1) ),
-	m_isoplanatic_angle_column( *this, boost::bind(&integrals::isoplanatic_angle,_1) ),
-	m_altitude_eff_column( *this, boost::bind(&integrals::altitude_eff,_1) ),
-	m_time_constant_column( *this, boost::bind(&integrals::time_constant,_1) ),
-	m_moment_0_column( *this, boost::bind(&integrals::moment_0,_1) ),
-	m_moment_2_column( *this, boost::bind(&integrals::moment_2,_1) ) {
-}
-
-bool integrals_regular_filter::operator() ( const integrals::integrals_sample& x ) const {
+	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 ),
+			    aints(0),
+			    aints(3) );
+}
+
+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_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 ) &&
@@ -212,6 +201,10 @@
 		is_regular( x.moment_0 ) &&
 		is_regular( x.moment_2 ) );
 }
-bool integrals_regular_filter::is_regular( double x ) const {

atmos-2.96.9.tar.gz/src/integrals.h -> atmos-2.97.1.tar.gz/src/integrals.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: integrals.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: integrals.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,86 +19,51 @@
 #ifndef _INTEGRALS_H
 #define _INTEGRALS_H
 
-#include <string>
-#include <stdexcept>
-
-#include <boost/tuple/tuple.hpp>
+#include <boost/function.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/vector.hpp>
+#include <boost/tuple/tuple.hpp>
 
-#include "table_storage.h"
+#include "indices.h"
 #include "nocase_less.h"
+#include "table_storage.h"
 
-#include "utils.h"
-#include "vector_generator.h"
-#include "weif.h"
-
+namespace atmos {
 namespace integrals {
-
-	struct integrals_sample {
-		double free_seeing;
-		double isoplanatic_angle;
-		double altitude_eff;
-		double time_constant;
-		double moment_0;
-		double moment_2;
+	using namespace boost::numeric::ublas;
+	
+	typedef double value_type;
+	typedef matrix< value_type > matrix_type;
+
+namespace detail {
+	struct sample {
+		typedef atmos::integrals::value_type value_type;
+
+		value_type free_seeing;
+		value_type isoplanatic_angle;
+		value_type altitude_eff;
+		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& time_constant,
+			const value_type& moment_0,
+			const value_type& moment_2 ):
+			free_seeing(free_seeing), isoplanatic_angle(isoplanatic_angle), altitude_eff(altitude_eff), time_constant(time_constant), moment_0(moment_0), moment_2(moment_2) {
+		}
 	};
+}
 
-	double free_seeing( const integrals_sample& x );
-	double isoplanatic_angle( const integrals_sample& x );
-	double altitude_eff( const integrals_sample& x );
-	double time_constant( const integrals_sample& x );
-	double moment_0( const integrals_sample& x );
-	double moment_2( const integrals_sample& x );
-};
-using namespace boost::numeric::ublas;
-
-class integrals_factory;
-class integrals_functor;
-
-class integrals_factory {
-	public:
-		typedef matrix< double > dcf_matrix_type;
-		typedef boost::function1< vector<double>, double > weight_type;
-		typedef std::map< std::string, dcf_matrix_type, utils::algo::nocase_less< std::string > > storage_type;
-		typedef integrals_functor functor_type;
-		typedef storage_type::const_iterator const_iterator;
-	private:
-		storage_type m_storage;
-
-		dcf_matrix_type calculate( const std::string& spectral_type, const weight_type& awf );
-	public:
-		integrals_factory();
-		integrals_factory( const integrals_factory& factory );
-		void swap( integrals_factory& factory ) throw();
-		integrals_factory& operator=( const integrals_factory& factory );
-		~integrals_factory();
-
-		functor_type get_functor( const std::string& spectral_type, const weight_type& awf );
-};
-
-class integrals_functor {
-	public:
-		typedef integrals_factory             factory_type;
-		typedef factory_type::dcf_matrix_type dcf_matrix_type;
-
-		typedef integrals::integrals_sample result_type;
-
-		typedef symmetric_matrix<double> indices_type;
-		typedef vector<double> indesis_type;
-	private:
-		const dcf_matrix_type& m_dcf;
-	public:
-		integrals_functor( const dcf_matrix_type& dcf ):
-			m_dcf( dcf ) {}
-
-		result_type operator() ( const indices_type& indices, const indesis_type& indesis, double mz ) const;
-};
-
-class integrals_storage:
-	public table_storage< integrals::integrals_sample > {
+class storage:
+	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 table_storage_column< integrals::integrals_sample, boost::function1< double, integrals::integrals_sample > > common_column_type;
 	typedef common_column_type free_seeing_column_type;
 	typedef common_column_type isoplanatic_angle_column_type;
 	typedef common_column_type altitude_eff_column_type;
@@ -113,7 +78,7 @@
 	moment_0_column_type          m_moment_0_column;
 	moment_2_column_type          m_moment_2_column;
 public:
-	integrals_storage( size_type reserve = 60 );
+	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; }
@@ -122,16 +87,61 @@
 	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; }
 
-	void push_back( const integrals::integrals_sample& integrals ) {
-		table_storage< integrals::integrals_sample >::push_back( integrals );
+	void push_back( const detail::sample& integrals ) {
+		table_storage< detail::sample >::push_back( integrals );
 	}
 };
 
-class integrals_regular_filter {
-public:
-	bool operator() ( const integrals::integrals_sample& x ) const;
-private:
-	bool is_regular( double x ) const;
+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;
+
+class functor_factory {
+	public:
+		typedef matrix_type dcf_matrix_type;
+		typedef boost::function1< vector<double>, double > weight_type;
+		typedef std::map< std::string, dcf_matrix_type, utils::algo::nocase_less< std::string > > storage_type;
+		typedef functor functor_type;
+	private:
+		storage_type m_storage;
+
+		dcf_matrix_type calculate( const std::string& spectral_type, const weight_type& awf );
+	public:
+		functor_factory();
+		functor_factory( const functor_factory& factory );
+		void swap( functor_factory& factory ) throw();
+		functor_factory& operator=( const functor_factory& factory );
+		~functor_factory();
+
+		functor_type operator() ( const std::string& spectral_type, const weight_type& awf );
 };
 
+class functor {
+	public:
+		typedef matrix_type dcf_matrix_type;
+
+		typedef atmos::indices::indices_type indices_type;
+		typedef atmos::indices::indesis_type indesis_type;
+		typedef atmos::integrals::detail::sample result_type;
+	private:
+		const dcf_matrix_type& m_dcf;
+	public:
+		functor( const dcf_matrix_type& dcf ):

atmos-2.96.9.tar.gz/src/io.cpp -> atmos-2.97.1.tar.gz/src/io.cpp Changed

@@ -32,7 +32,7 @@
 		stm.iword(iword_index) = true;
 		return stm;
 	}
-	std::ostream& operator<<( std::ostream& stm, const datum& m ){
+	std::ostream& operator<<( std::ostream& stm, const atmos::datum& m ){
 		boost::io::ios_flags_saver fsaver( stm );
 		if( std::ostream::sentry( stm ) ){
 			int w = ( stm.flags() & std::ios::scientific )? 8 : 5;
@@ -43,13 +43,13 @@
 		}
 		return stm;
 	}
-	std::ostream& operator<<( std::ostream& stm, const object& obj ){
+	std::ostream& operator<<( std::ostream& stm, const atmos::object& obj ){
 		if( std::ostream::sentry( stm ) ){
 			stm << obj.hd << stm.fill() << obj.name << stm.fill() << std::hours << std::sexagesimal << yalpa::make_angle( obj.coordinates.ra() ) << std::degrees << stm.fill() << yalpa::make_angle( obj.coordinates.dec() ) << stm.fill() << obj.vmag << stm.fill() << obj.bvind << stm.fill() << obj.sclass;
 		}
 		return stm;
 	}
-	std::istream& operator>>( std::istream& stm, object& obj ){
+	std::istream& operator>>( std::istream& stm, atmos::object& obj ){
 		if( std::istream::sentry( stm ) ){
 			//O 2009-07-19 18:04:29 5191 Eta_UMa 13 47 32 +49 18 48  1.86 -0.19 B55
 			double ra = 0,dec = 0;

atmos-2.96.9.tar.gz/src/io.h -> atmos-2.97.1.tar.gz/src/io.h Changed

@@ -50,23 +50,29 @@
 	template<class T> std::ostream& operator<<( std::ostream& stm, const std::vector<T>& x ) {
 		typedef typename std::vector<T>::const_iterator const_iterator;
 
-		for( const_iterator it = x.begin(); it != x.end(); ++it ){
-			stm << *it << stm.fill();
+		const_iterator it = x.begin();
+		if( it != x.end() )
+			stm << *(it++);
+		for( ; it != x.end(); ++it ){
+			stm << stm.fill() << *it;
 		}
 		return stm;
 	}
 	template<class T, class VA> std::ostream& operator<<( std::ostream& stm, const boost::numeric::ublas::vector<T,VA>& x ) {
 		typedef typename boost::numeric::ublas::vector<T,VA>::const_iterator const_iterator;
 
-		for( const_iterator it = x.begin(); it != x.end(); ++it ){
-			stm << *it << stm.fill();
+		const_iterator it = x.begin();
+		if( it != x.end() )
+			stm << *(it++);
+		for( ; it != x.end(); ++it ){
+			stm << stm.fill() << *it;
 		}
 		return stm;
 	}
 	std::ostream& indic( std::ostream& stm );
-	std::ostream& operator<<( std::ostream& stm, const datum& m );
-	std::ostream& operator<<( std::ostream& stm, const object& obj );
-	std::istream& operator>>( std::istream& stm, object& obj );
+	std::ostream& operator<<( std::ostream& stm, const atmos::datum& m );
+	std::ostream& operator<<( std::ostream& stm, const atmos::object& obj );
+	std::istream& operator>>( std::istream& stm, atmos::object& obj );
 }
 
 #endif // _IO_H

atmos-2.97.1.tar.gz/src/local_context.h Added

@@ -0,0 +1,41 @@
+/*
+    $Id: local_context.h 180 2011-01-30 15:56:33Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _LOCAL_CONTEXT_H
+#define _LOCAL_CONTEXT_H
+
+namespace utils {
+
+template<class Context> class local_context {
+	public:
+		typedef Context                context_type;
+		typedef context_type&          context_reference_type;
+		typedef const context_type&    context_const_reference_type;
+		typedef local_context<Context> local_context_type;
+	private:
+		context_reference_type m_ctx;
+	protected:
+		inline context_reference_type       context()       throw() { return m_ctx; }
+		inline context_const_reference_type context() const throw() { return m_ctx; }
+	public:
+		local_context( context_reference_type ctx ) throw() : m_ctx(ctx) {}
+};
+
+}
+
+#endif // _LOCAL_CONTEXT_H

atmos-2.96.9.tar.gz/src/main.cpp -> atmos-2.97.1.tar.gz/src/main.cpp Changed

atmos-2.97.1.tar.gz/src/params.cpp Added

@@ -0,0 +1,132 @@
+/*
+    $Id: params.cpp 185 2011-01-31 12:01:43Z matwey $
+    Copyright (C) 2009  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <stdexcept>
+
+#include <sstream>
+
+#include <yalpa/spheric.h>
+#include <yalpa/io.h>
+
+#include "atmos.h"
+#include "params.h"
+
+namespace atmos {
+namespace parameters {
+
+container::container() {}
+container::~container() {}
+bool container::key_equal( const key_type& k1, const key_type& k2 ) {
+	return !m_storage.key_comp()( k1, k2 ) && !m_storage.key_comp()( k2, k1 );
+}
+bool container::value_equal( const value_type& v1, const value_type& v2 ){
+	return !m_storage.value_comp()( v1, v2 ) && !m_storage.value_comp()( v2, v1 );
+}
+bool container::has(const storage_type::key_type& key) const {
+	storage_type::const_iterator it = m_storage.find( key );
+	return ( it != m_storage.end() );
+}
+const container::storage_type::mapped_type& container::at(const storage_type::key_type& key) const {
+	storage_type::const_iterator it = m_storage.find( key );
+	if( it != m_storage.end() )
+		return it->second;
+	throw std::range_error( std::string("Not found value for key \"") + key + std::string("\""));
+}
+
+proxy::proxy( const container& container ): m_container( container ) {}
+proxy::~proxy() {}
+
+value_type proxy::longitude() const {
+	value_type longitude = 0.0;
+	std::istringstream iss( m_container.at("general/site/longitude") );
+	iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
+	iss >> std::sexagesimal >> std::hours >> yalpa::make_angle( longitude );
+	return longitude;
+}
+value_type proxy::latitude() const {
+	value_type latitude = 0.0;
+	std::istringstream iss( m_container.at("general/site/latitude") );
+	iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
+	iss >> std::sexagesimal >> std::degrees >> yalpa::make_angle( latitude );
+	return latitude;
+}
+yalpa::geographic< value_type > proxy::site() const {
+	return yalpa::geographic< value_type >( longitude(), latitude() );
+}
+
+vector_type proxy::inner_diam()    const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as<value_type>("optics/segment a/inner");
+	ret(1) = m_container.as<value_type>("optics/segment b/inner");
+	ret(2) = m_container.as<value_type>("optics/segment c/inner");
+	ret(3) = m_container.as<value_type>("optics/segment d/inner");
+	return ret * ( magnification()*0.1 );
+}
+vector_type proxy::outer_diam()    const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as<value_type>("optics/segment a/outer");
+	ret(1) = m_container.as<value_type>("optics/segment b/outer");
+	ret(2) = m_container.as<value_type>("optics/segment c/outer");
+	ret(3) = m_container.as<value_type>("optics/segment d/outer");
+	return ret * ( magnification()*0.1 );
+}
+vector_type proxy::nonlinearity() const {
+	typedef vector_type::value_type value_type;
+	std::vector< std::string > ctr( counters() );
+	vector_type ret( ctr.size() );
+	for( std::vector< std::string >::size_type i = 0; i < ctr.size(); ++i ) {
+		ret(i) = m_container.as< value_type >( (boost::format("Devices/Bicounter %1%/NonLinearity%2%") % ctri0 % ctri1).str() );
+	}
+	return ret * 1e-6 / exposition();
+}
+vector_type proxy::nonpoisson() const {
+	typedef vector_type::value_type value_type;
+	std::vector< std::string > ctr( counters() );
+	vector_type ret( ctr.size() );
+	for( std::vector< std::string >::size_type i = 0; i < ctr.size(); ++i ) {
+		ret(i) = m_container.as< value_type >( (boost::format("Devices/Bicounter %1%/NonPoisson%2%") % ctri0 % ctri1).str() );
+	}
+	return ret;
+}
+vector_type proxy::scatter() const {
+	typedef vector_type::value_type value_type;
+	vector_type ret(atmos::channels);
+	ret(0) = m_container.as< value_type >("channels/scatter/channela",0.0);
+	ret(1) = m_container.as< value_type >("channels/scatter/channelb",0.0);
+	ret(2) = m_container.as< value_type >("channels/scatter/channelc",0.0);
+	ret(3) = m_container.as< value_type >("channels/scatter/channeld",0.0);
+	return ret;
+}
+std::vector< std::string > proxy::counters() const {
+	std::vector< std::string > ret;
+	ret.resize(atmos::channels);
+	ret0 = m_container.at("channels/configuration/channela");
+	ret1 = m_container.at("channels/configuration/channelb");
+	ret2 = m_container.at("channels/configuration/channelc");
+	ret3 = m_container.at("channels/configuration/channeld");
+	for( std::vector< std::string >::iterator it = ret.begin(); it != ret.end(); ++it ){
+		//FIXME: use regex
+		if( it->size() < 10 || it->substr(0,7) != "Counter" || ( (*it)8 != '1' && (*it)8 != '2' ) || ( (*it)9 != 'A' && (*it)9 != 'B' ) ) throw std::logic_error( std::string("Bad value in channel configuration: ") + *it);
+		*it = it->substr(8,2);
+	}
+	return ret;
+}
+	
+}}

atmos-2.96.9.tar.gz/src/params.h -> atmos-2.97.1.tar.gz/src/params.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: params.h 98 2010-03-01 12:55:46Z matwey $
+    $Id: params.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,61 +19,44 @@
 #ifndef _PARAMS_H
 #define _PARAMS_H
 
-#include <algorithm>
-#include <cmath>
+#include <map>
 #include <string>
-#include <exception>
-#include <stdexcept>
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/lexical_cast.hpp>
 #include <boost/format.hpp>
+#include <boost/lexical_cast.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
-#include <yalpa/spheric.h>
-#include <yalpa/io.h>
+#include "nocase_less.h"
 
-using namespace boost::numeric;
-// use ublas::vector for ublas vector and std::vector for std vector
+namespace atmos {
+namespace parameters {
 
-/**
-  @brief Container for handling parameters extracted from MASS file.
-*/
-template<class ContainerType> class parameters_container {
+using namespace boost::numeric::ublas;
+
+class container {
 public:
-	typedef ContainerType storage_type;
-	typedef typename storage_type::key_type   key_type;
-	typedef typename storage_type::value_type value_type;
-	typedef typename storage_type::const_iterator const_iterator;
+	typedef std::map< std::string, std::string, utils::algo::nocase_less< std::string > > storage_type;
+	typedef storage_type::key_type   key_type;
+	typedef storage_type::value_type value_type;
+	typedef storage_type::const_iterator const_iterator;
 private:
 	storage_type m_storage;
 public:
-	parameters_container() {}
-	~parameters_container() {}
+	container();
+	~container();
 
-	bool key_equal( const key_type& k1, const key_type& k2 ) {
-		return !m_storage.key_comp()( k1, k2 ) && !m_storage.key_comp()( k2, k1 );
-	}
-	bool value_equal( const value_type& v1, const value_type& v2 ){
-		return !m_storage.value_comp()( v1, v2 ) && !m_storage.value_comp()( v2, v1 ) ;
-	}
-
-	bool has(const typename storage_type::key_type& key) const {
-		typename storage_type::const_iterator it = m_storage.find( key );
-		return ( it != m_storage.end() );
-	}
-	const typename storage_type::mapped_type& at(const typename storage_type::key_type& key) const {
-		typename storage_type::const_iterator it = m_storage.find( key );
-		if( it != m_storage.end() )
-			return it->second;
-		throw std::range_error( std::string("Not found value for key \"") + key + std::string("\""));
-	}
-	typename storage_type::mapped_type& operator (const typename storage_type::key_type& key) {
+	bool key_equal( const key_type& k1, const key_type& k2 );
+	bool value_equal( const value_type& v1, const value_type& v2 );
+	bool has(const storage_type::key_type& key) const;
+	const storage_type::mapped_type& at(const storage_type::key_type& key) const;
+	
+	storage_type::mapped_type& operator (const storage_type::key_type& key) {
 		return m_storagekey;
 	}
-	template<class Y> Y as(const typename storage_type::key_type& key) const {
+	template<class Y> Y as(const storage_type::key_type& key) const {
 		return boost::lexical_cast<Y>( at(key) );
 	}
-	template<class Y> Y as(const typename storage_type::key_type& key, const Y& def ) const {
+	template<class Y> Y as(const storage_type::key_type& key, const Y& def ) const {
 		if( has( key ) )
 			return boost::lexical_cast<Y>( at(key) );
 		return def;
@@ -81,110 +64,43 @@
 
 	const_iterator begin() const { return m_storage.begin(); }
 	const_iterator end() const { return m_storage.end(); }
-
 };
 
+typedef double               value_type;
+typedef vector< value_type > vector_type;
 
-/**
-  @brief Proxy for access to parameters extracted from MASS file and stored in container.
-*/
-template<class ParamType> class parameters_proxy {
-public:
-	typedef ParamType params_type;
+class proxy {
 private:
-	const params_type& m_container;
+	const container& m_container;
 public:
-	parameters_proxy( const params_type& container ): m_container( container ) {}
-	~parameters_proxy() {}
+	proxy( const container& container );
+	~proxy();
 
 	std::string sitename() const { return m_container.at("general/site/sitename"); }
-	double longitude() const {
-		double longitude = 0.0;
-		std::istringstream iss( m_container.at("general/site/longitude") );
-		iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
-		iss >> std::sexagesimal >> std::hours >> yalpa::make_angle( longitude );
-		return longitude;
-	}
-	double latitude() const {
-		double latitude = 0.0;
-		std::istringstream iss( m_container.at("general/site/latitude") );
-		iss.exceptions(std::ios_base::badbit | std::ios_base::failbit);
-		iss >> std::sexagesimal >> std::degrees >> yalpa::make_angle( latitude );
-		return latitude;
-	}
-	yalpa::geographic<double> site() const {
-		return yalpa::geographic<double>( longitude(), latitude() );
-	}
-
-	double accumtime()     const { return m_container.template as<double>("operations/normal mode/accumtime"); }
-	double basetime()      const { return m_container.template as<double>("operations/normal mode/basetime"); }
-	double exposition()    const { return m_container.template as<double>("operations/normal mode/exposition"); }
-
-	double magnification() const { return m_container.template as<double>("optics/common/magnification"); }
-
-	//FIXME: make function with wildcards
-	ublas::vector<double> inner_diam()    const {
-		ublas::vector<double> ret(4);
-		ret(0) = m_container.template as<double>("optics/segment a/inner");
-		ret(1) = m_container.template as<double>("optics/segment b/inner");
-		ret(2) = m_container.template as<double>("optics/segment c/inner");
-		ret(3) = m_container.template as<double>("optics/segment d/inner");
-		return ret * ( magnification()*0.1 );
-	}
-	ublas::vector<double> outer_diam()    const {
-		ublas::vector<double> ret(4);
-		ret(0) = m_container.template as<double>("optics/segment a/outer");
-		ret(1) = m_container.template as<double>("optics/segment b/outer");
-		ret(2) = m_container.template as<double>("optics/segment c/outer");
-		ret(3) = m_container.template as<double>("optics/segment d/outer");
-		return ret * ( magnification()*0.1 );
-	}
-	ublas::vector<double> nonlinearity() const {
-		std::vector< std::string > ctr( counters() );
-		ublas::vector<double> ret( ctr.size() );
-		for( typename std::vector< std::string >::const_iterator it = ctr.begin(); it != ctr.end(); ++it ) {
-			ret( it - ctr.begin() ) = m_container.template as<double>( (boost::format("Devices/Bicounter %1%/NonLinearity%2%") % (*it)0 % (*it)1).str() );
-		}
-		return ret * 1e-6 / exposition();
-	}
-	ublas::vector<double> nonpoisson() const {
-		std::vector< std::string > ctr( counters() );
-		ublas::vector<double> ret( ctr.size() );
-		for( typename std::vector< std::string >::const_iterator it = ctr.begin(); it != ctr.end(); ++it ) {
-			ret( it - ctr.begin() ) = m_container.template as<double>( (boost::format("Devices/Bicounter %1%/NonPoisson%2%") % (*it)0 % (*it)1).str() );
-		}
-		return ret;
-	}
-	ublas::vector<double> scatter() const {
-		ublas::vector<double> ret( 4 );
-		ret(0) = m_container.template as<double>("channels/scatter/channela",0.0);
-		ret(1) = m_container.template as<double>("channels/scatter/channelb",0.0);
-		ret(2) = m_container.template as<double>("channels/scatter/channelc",0.0);
-		ret(3) = m_container.template as<double>("channels/scatter/channeld",0.0);
-		return ret;
-	}
-	//FIXME: use tuple
-	std::vector< std::string > counters() const {
-		std::vector< std::string > ret;
-		ret.resize(4);
-		ret0 = m_container.at("channels/configuration/channela");
-		ret1 = m_container.at("channels/configuration/channelb");
-		ret2 = m_container.at("channels/configuration/channelc");
-		ret3 = m_container.at("channels/configuration/channeld");
-		for( typename std::vector< std::string >::iterator it = ret.begin(); it != ret.end(); ++it ){
-			//FIXME: use regex
-			if( it->size() < 10 || it->substr(0,7) != "Counter" || ( (*it)8 != '1' && (*it)8 != '2' ) || ( (*it)9 != 'A' && (*it)9 != 'B' ) ) throw std::logic_error( std::string("Bad value in channel configuration: ") + *it);
-			*it = it->substr(8,2);
-		}
-		return ret;
-	}
-

atmos-2.96.9.tar.gz/src/photometry.h -> atmos-2.97.1.tar.gz/src/photometry.h Changed

atmos-2.96.9.tar.gz/src/profile.cpp -> atmos-2.97.1.tar.gz/src/profile.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: profile.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: profile.cpp 182 2011-01-31 11:26:03Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,101 +16,137 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <lsp/nnls.h>
-#include <boost/bind.hpp>
+#include <algorithm>
+
 #include <boost/format.hpp>
+#include <boost/mem_fn.hpp>
+#include <boost/numeric/ublas/banded.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
-#include <boost/numeric/ublas/io.hpp>
-#include <iostream>
-#include <iomanip>
-#include <algorithm>
-#include <numeric>
 
-#include "io.h"
+#include <lsp/cholesky_decomposition.h>
+#include <lsp/nnls.h>
+
+#include "atmos.h"
+#include "apply_to_all.h"
+#include "covariance_decomposition.h"
 #include "profile.h"
 #include "utils.h"
 
-const double MAXCHI2 = 999.0;
-const double JPOW_SEE = 0.6 ;
-const double LPOW_SEE = -0.2;
-const double KS = 1.73588e+07 ;
-
-using namespace boost::numeric::ublas;
+namespace atmos {
+namespace mixture {
 
-restoration_functor::restoration_functor( const vector<double>& grid, const symmetric_matrix< double >& corr_matrix, double secz, const weight_type& basewf ):
-	m_grid( grid ), m_cmatrix(corr_matrix), m_secz(secz), m_lambda_eff( 0.5 ) {
-	using namespace utils::ublas;
+namespace {
+	
+vector_type reorder_flat( const atmos::indices::indices_type& indices ) {
+	assert( indices.size1() == atmos::channels );
+	assert( indices.size2() == atmos::channels );
+	
+	vector_type inds( atmos::indices_size );
+	inds(0) = indices(0,0); inds(1) = indices(1,1); inds(2) = indices(2,2); inds(3) = indices(3,3);
+	inds(4) = indices(0,1); inds(5) = indices(0,2); inds(6) = indices(0,3);
+	inds(7) = indices(1,2); inds(8) = indices(1,3); inds(9) = indices(2,3);
+	return inds;
+}
 
-	if( zero_elements( diag( m_cmatrix ) ) > 0 ) {
-		throw std::logic_error( (boost::format("Zero variances: %1% x %2% %3%") % m_cmatrix.size1() % m_cmatrix.size2()  % vector<double>( diag( m_cmatrix )) ).str().c_str() );
+}
+	
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_mixture_column( *this, boost::mem_fn(&detail::sample::mixture) ) {
+}
+storage::storage( const atmos::indices::storage& indices, const atmos::dimm::storage& dimm, double dimmscale ):
+	table_storage< detail::sample >( indices.size() ),
+	m_mixture_column( *this, boost::mem_fn(&detail::sample::mixture) ) {
+	using namespace boost::numeric::ublas;
+	
+	if( dimm.size() > 1 ) {
+		double dimmscale2 = dimmscale*dimmscale;
+		dimm::storage::const_iterator dimm_current = dimm.begin();
+		dimm::storage::const_iterator dimm_next = dimm_current + 1;
+		indices::storage::const_iterator it = indices.begin();
+
+		while( dimm_current != dimm.end() ) {
+			while( it != indices.end() ) {
+				if( dimm_next != dimm.end() && dimm_next->timestamp <= it->timestamp )
+					break;
+
+				mixture_type mixture( atmos::indices_size + atmos::dimm_indices_size );
+				project( mixture, range(0,atmos::indices_size) ) = reorder_flat(it->indices);
+				mixture(10) = dimmscale2 * dimm_current->xvars(1) * 1e+11;
+				mixture(11) = dimmscale2 * dimm_current->yvars(1) * 1e+11;
+				++it;
+				push_back( mixture );
+			}
+			dimm_current = dimm_next;
+			++dimm_next;
+		}
+	} else {
+		for( indices::storage::const_iterator it = indices.begin(); it != indices.end(); ++it ) {
+			push_back( reorder_flat(it->indices) );
+		}
 	}
+}
 
-	lsp::cholesky_decomposition( m_cmatrix, upper_tag() );
-	triangular_adaptor< symmetric_matrix< double >, lower > R( m_cmatrix );
-	m_wmatrix = identity_matrix< double >( corr_matrix.size1() );
-	inplace_solve( R, m_wmatrix, lower_tag() );
+}
+namespace profile {
 
-	if( zero_elements( diag( m_wmatrix ) ) > 0 ) {
-		throw std::logic_error( (boost::format("Zero weights: %1%") % vector<double>( diag( m_wmatrix )) ).str().c_str() );
-	}
+const double MAXCHI2 = 999.0;
+const double JPOW_SEE = 0.6 ;
+const double LPOW_SEE = -0.2;
+const double KL = 1.1486983549970351; // 0.5 ** LPOW_SEE
+const double KS = 1.73588e+07;
 
-	m_pmatrix = prod( m_wmatrix, subrange( vector_of_vector_to_matrix( project(grid*secz, basewf) ) ,0,m_wmatrix.size2(),0,grid.size()) );
-}
-vector< double > restoration_functor::flat_reorder( const matrix<double>& indices ) {
-	vector< double > inds(10);
-	inds(0) = indices(0,0); inds(1) = indices(1,1); inds(2) = indices(2,2); inds(3) = indices(3,3);
-	inds(4) = indices(0,1); inds(5) = indices(0,2); inds(6) = indices(0,3);
-	inds(7) = indices(1,2); inds(8) = indices(1,3); inds(9) = indices(2,3);
-	return inds;
+functor::functor( const grid_type& grid, const mixture_covariance_type& cov_matrix, double secz, const weight_type& basewf ):
+	m_grid( grid ),
+	m_weight_matrix( utils::ublas::aetkin_weight_matrix(cov_matrix) ),
+	m_invweight_matrix( solve( m_weight_matrix, identity_matrix< value_type >( m_weight_matrix.size1() ), lower_tag() ) ),
+	m_secz( secz ) {
+	using namespace utils::ublas;
+	using namespace std;
+
+	m_problem_matrix = prod( m_weight_matrix, subrange( vector_of_vector_to_matrix( project(grid*secz, basewf) ), 0, m_weight_matrix.size2(),0,grid.size()) );
 }
-restoration_functor::result_type restoration_functor::operator() ( const vector<double>& indices ) const {
-	vector< double > inds = prod( m_wmatrix, indices );
+
+functor::result_type functor::operator() ( const mixture_type& indices ) const {
+	vector_type weighted_indices = prod( m_weight_matrix, indices );
 
 	result_type ret;
 
-	matrix< double > cov = zero_matrix<double>( m_pmatrix.size2(), m_pmatrix.size2() );
+	matrix_type cov = zero_matrix<double>( m_problem_matrix.size2(), m_problem_matrix.size2() );
 
-	lsp::nnls< matrix< double >, vector< double > > nnls( m_pmatrix, inds );
+	lsp::nnls< matrix_type, vector_type > nnls( m_problem_matrix, weighted_indices );
 	nnls.solve( ret.cn2, cov );
 
-	ret.residual = prod( m_pmatrix, ret.cn2 ) - inds;
+	ret.residual = prod( m_problem_matrix, ret.cn2 ) - weighted_indices;
 
 	ret.chi2 = std::min( inner_prod( ret.residual, ret.residual ), MAXCHI2 );
 
 	ret.cn2 /= m_secz;
 	return ret;
 }
-double restoration_functor::seeing( const vector<double>& cn2 ) const {
-	double cn2_sum = sum( cn2 );
-	double secz = 1.0; // already corrected cn2 at input
-	return KS*pow( cn2_sum/secz, JPOW_SEE )*pow( m_lambda_eff, LPOW_SEE );
+double functor::seeing( const cn2_type& cn2 ) const {
+	return KS*KL*std::pow( sum( cn2 ), JPOW_SEE );
 }
-vector<double> restoration_functor::indices_residual( const vector<double>& residuals ) const {
-	return prod( triangular_adaptor< const matrix< double >, lower >(m_cmatrix), residuals );
+vector_type functor::indices_residual( const residual_type& residual ) const {
+	return prod( m_invweight_matrix, residual );
 }
 
-profile::chi2_type profile::chi2( const profile_sample& x ) {
-	return x.chi2;
-}
-const profile::cn2_type& profile::cn2( const profile_sample& x ) {
-	return x.cn2;
-}
-const profile::residual_type& profile::residual( const profile_sample& x ) {
-	return x.residual;
+storage::storage( size_type reserve ):
+	table_storage< detail::sample >( reserve ),
+	m_cn2_column( *this, boost::mem_fn(&detail::sample::cn2) ),
+	m_residual_column( *this, boost::mem_fn(&detail::sample::residual)  ),
+	m_chi2_column( *this, boost::mem_fn(&detail::sample::chi2) ) {
 }
 
-profile_storage::profile_storage( size_type reserve ):
-	table_storage< profile::profile_sample >( reserve ),
-	m_cn2_column( *this, boost::bind(&profile::cn2, _1) ),
-	m_residual_column( *this, boost::bind(&profile::residual, _1)  ),
-	m_chi2_column( *this, boost::bind(&profile::chi2, _1) )
-{
-}
+namespace filter {
 
-chi2_gt_than_filter::chi2_gt_than_filter( double chi2 ): chi2_(chi2) {
+chi2_gt_than::chi2_gt_than( const value_type& chi2 ): chi2_(chi2) {
 }

atmos-2.96.9.tar.gz/src/profile.h -> atmos-2.97.1.tar.gz/src/profile.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: profile.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: profile.h 182 2011-01-31 11:26:03Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,101 +19,146 @@
 #ifndef _PROFILE_H
 #define _PROFILE_H
 
-#include <string>
-#include <stdexcept>
-
-#include <lsp/cholesky_decomposition.h>
-
+#include <boost/function.hpp>
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/vector.hpp>
 #include <boost/numeric/ublas/symmetric.hpp>
 
+#include "indices.h"
+#include "dimm.h"
 #include "table_storage.h"
-#include "weif.h"
 
-namespace profile {
+namespace atmos {
+namespace mixture {
 	using namespace boost::numeric::ublas;
 
 	typedef double                 value_type;
 	typedef vector<value_type>     vector_type;
 	typedef vector_type::size_type size_type;
+	typedef symmetric_matrix<value_type> symmetric_matrix_type;
 
-	typedef vector_type cn2_type;
-	typedef vector_type residual_type;
-	typedef value_type  chi2_type;
+	typedef vector_type           mixture_type;
+	typedef symmetric_matrix_type mixture_covariance_type;
 
-	struct profile_sample {
-		cn2_type       cn2;
-		residual_type  residual;
-		chi2_type      chi2;
-	};
+namespace detail {
+	struct sample {
+		typedef atmos::mixture::mixture_type mixture_type;
 
-	chi2_type chi2( const profile_sample& x );
-	const cn2_type& cn2( const profile_sample& x );
-	const residual_type& residual( const profile_sample& x );
-};
+		mixture_type mixture;
 
-template<class VE> typename VE::size_type zero_elements( const vector_expression<VE>&x ) {
-    typename VE::size_type n = 0;
-    for( typename VE::size_type j = 0; j < x().size(); j++ )
-	    n += ( x()(j) == 0 || std::isinf(x()(j)) || std::isnan(x()(j)) ) ? 1 : 0;
-    return n;
+		sample( const mixture_type& mixture ): mixture(mixture) {
+		}
+	};
 }
 
-class restoration_functor {
-	public:
-		typedef profile::profile_sample result_type;
+class storage:
+	public table_storage< detail::sample > {
+public:
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::mixture_type > mixture_column_type;
+private:
+	mixture_column_type m_mixture_column;
+public:
+	storage( size_type reserve = 60 );
+	storage( const atmos::indices::storage& indices, const atmos::dimm::storage& dimm, double dimmscale );
 
-		typedef vector<double> grid_type;
-		typedef matrix<double> weight_matrix;
-		typedef boost::function1< vector<double>, double > weight_type;
-	private:
-		grid_type     m_grid;
-		weight_matrix m_wmatrix;
-		symmetric_matrix< double > m_cmatrix;
-		matrix< double > m_pmatrix;
-		double        m_secz;
-		double        m_lambda_eff;
+	const mixture_column_type& mixture() const { return m_mixture_column; }
 
-	public:
-		static vector< double > flat_reorder( const matrix<double>& indices );
+	void push_back( const mixture_type& mixture ) {
+		table_storage< detail::sample >::push_back( detail::sample( mixture ) );
+	}
+};
+}
+namespace profile {
+	using namespace boost::numeric::ublas;
 
-		restoration_functor( const vector<double>& grid, const symmetric_matrix< double >& corr_matrix, double secz, const weight_type& basewf );
+	typedef double                 value_type;
+	typedef vector<value_type>     vector_type;
+	typedef matrix<value_type>     matrix_type;
+	typedef triangular_matrix<value_type> triangular_matrix_type;
+	typedef vector_type::size_type size_type;
 
-		result_type operator() ( const vector<double>& indices ) const;
+	typedef vector_type cn2_type;
+	typedef vector_type residual_type;
+	typedef value_type  chi2_type;
 
-		double seeing( const vector<double>& cn2 ) const;
-		vector<double> indices_residual( const vector<double>& residuals ) const;
-};
+namespace detail {
+	struct sample {
+		typedef atmos::profile::cn2_type      cn2_type;
+		typedef atmos::profile::residual_type residual_type;
+		typedef atmos::profile::chi2_type     chi2_type;
+	
+		cn2_type       cn2;
+		residual_type  residual;
+		chi2_type      chi2;
+
+		sample( const cn2_type& cn2, const residual_type& residual, const chi2_type& chi2 ):
+			cn2(cn2), residual(residual), chi2(chi2) {
+		}
+		sample() {
+		}
+	};
+}
 
-class profile_storage:
-	public table_storage< profile::profile_sample > {
+class storage:
+	public table_storage< detail::sample > {
 public:
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::cn2_type,      profile::profile_sample > > cn2_column_type;
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::residual_type, profile::profile_sample > > residual_column_type;
-	typedef table_storage_column< profile::profile_sample, boost::function1< profile::chi2_type,     profile::profile_sample > > chi2_column_type;
+	typedef storage self_type;
+	typedef table_storage_column_getter< self_type, detail::sample::cn2_type >      cn2_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::residual_type > residual_column_type;
+	typedef table_storage_column_getter< self_type, detail::sample::chi2_type >     chi2_column_type;
 private:
 	cn2_column_type      m_cn2_column;
 	residual_column_type m_residual_column;
 	chi2_column_type     m_chi2_column;
 public:
-	profile_storage( size_type reserve = 60 );
+	storage( size_type reserve = 60 );
 
 	const cn2_column_type&      cn2()      const { return m_cn2_column; }
 	const residual_column_type& residual() const { return m_residual_column; }
 	const chi2_column_type&     chi2()     const { return m_chi2_column; }
 
-	void push_back( const profile::profile_sample& profile ) {
-		table_storage< profile::profile_sample >::push_back( profile );
+	void push_back( const detail::sample& profile ) {
+		table_storage< detail::sample >::push_back( profile );
 	}
 };
 
-class chi2_gt_than_filter {
+class functor {
+	public:
+		typedef detail::sample result_type;
+		typedef atmos::mixture::mixture_type mixture_type;
+		typedef atmos::mixture::mixture_covariance_type mixture_covariance_type;
+
+		typedef vector_type grid_type;
+		typedef boost::function1< vector<double>, double > weight_type;
 	private:
-		double chi2_;
+		grid_type     m_grid;
+		triangular_matrix_type m_weight_matrix;
+		triangular_matrix_type m_invweight_matrix;
+		matrix_type   m_problem_matrix;
+		double        m_secz;
+
 	public:
-		chi2_gt_than_filter( double chi2 );
-		bool operator()( const profile::profile_sample& x ) const;
+		functor( const grid_type& grid, const mixture_covariance_type& cov_matrix, double secz, const weight_type& basewf );
+
+		result_type operator() ( const mixture_type& indices ) const;
+
+		double seeing( const cn2_type& cn2 ) const;
+		vector_type indices_residual( const residual_type& residual ) const;
 };
 
+namespace filter {
+	class chi2_gt_than {
+	public:
+		typedef atmos::profile::detail::sample sample_type;
+		typedef atmos::profile::value_type value_type;

atmos-2.96.9.tar.gz/src/spectral_factory.cpp -> atmos-2.97.1.tar.gz/src/spectral_factory.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: spectral_factory.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: spectral_factory.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,16 +22,15 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/filesystem/operations.hpp>
 
-#include <boost/numeric/ublas/io.hpp>
-
+#include "io.h"
+#include "spectral_factory.h"
 #include "weif.h"
 #include "utils.h"
 
-#include "spectral_factory.h"
-
-using namespace boost::numeric::ublas;
+namespace atmos {
+namespace weight {
 
-std::pair< spectral::vector_type, spectral::vector_type > spectral::load_reaction( const boost::filesystem::path& file ) {
+std::pair< vector_type, vector_type > load_reaction( const boost::filesystem::path& file ) {
 	std::ifstream fstm( file.string().c_str() );
 	if( fstm.fail() )
 		throw std::runtime_error( "Can't open file " + file.string() );
@@ -79,7 +78,7 @@
 	return ret;
 }
 
-std::pair< spectral::vector_type, spectral::vector_type > spectral::photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response ) {
+std::pair< vector_type, vector_type > photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response ) {
 	std::pair< vector_type, vector_type > ret;
 
 	assert( spectral_lambda.size() == spectral_response.size() );
@@ -112,7 +111,7 @@
 	return ret;
 }
 
-void spectral::dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid ) {
+void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid ) {
 	static const char* name = { "A", "B", "C", "D", "AB", "AC", "AD", "BC", "BD", "CD", "L", "T" };
 	static const unsigned int name_size = sizeof(name) / sizeof(name0);
 	const double wscale = 1e+11;
@@ -146,37 +145,38 @@
 	}
 }
 
-spectral_factory::spectral_factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response ):
+factory::factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response ):
 	m_mass_response(mass_response),
 	m_ccd_response(ccd_response),
 	m_wdir(wd),
 	m_basis_grid( lg_generator(0.0, 0.04, 0.25), 50 ),
-	m_mass_reaction( spectral::load_reaction(mass_response) ),
-	m_dimm_reaction( spectral::load_reaction(ccd_response) ) {
+	m_mass_reaction( load_reaction(mass_response) ),
+	m_dimm_reaction( load_reaction(ccd_response) ) {
 
 	if( ! (boost::filesystem::exists( wd ) && boost::filesystem::is_directory( wd )) )
 		throw std::runtime_error( wd.string() + " is not a directory or doesn't exists");
 }
-spectral_factory::~spectral_factory() {
-
+factory::~factory() {
 }
-const boost::filesystem::path spectral_factory::get_spectrum_filename( const std::string& spectype ) const {
+const boost::filesystem::path factory::get_spectrum_filename( const std::string& spectype ) const {
 	std::string spectral_type( spectype );
 	boost::to_lower( spectral_type );
 	boost::filesystem::path spect( m_wdir / (spectral_type + ".sp") );
 	return spect;
 }
-const spectral_factory::weight_type& spectral_factory::spectral( const std::string& spectype, const vector<double>& inner_diam, const vector<double>& outer_diam, double dimmbase, double dimmsize ) {
+const factory::weight_type& factory::weight( const std::string& spectype, const vector<double>& inner_diam, const vector<double>& outer_diam, double dimmbase, double dimmsize ) {
 	storage_type::const_iterator it = m_storage.find( spectype );
 	if( it != m_storage.end() )
 		return it->second;
 
-	std::pair< vector_type, vector_type > spectral_response = spectral::load_reaction( get_spectrum_filename( spectype ) );
+	std::pair< vector_type, vector_type > spectral_response = load_reaction( get_spectrum_filename( spectype ) );
 
-	std::pair< vector_type, vector_type > mass_response = spectral::photon_distribution( spectral_response.first, spectral_response.second, m_mass_reaction.first, m_mass_reaction.second );
-	std::pair< vector_type, vector_type > dimm_response = spectral::photon_distribution( spectral_response.first, spectral_response.second, m_dimm_reaction.first, m_dimm_reaction.second );
+	std::pair< vector_type, vector_type > mass_response = photon_distribution( spectral_response.first, spectral_response.second, m_mass_reaction.first, m_mass_reaction.second );
+	std::pair< vector_type, vector_type > dimm_response = photon_distribution( spectral_response.first, spectral_response.second, m_dimm_reaction.first, m_dimm_reaction.second );
 
-	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectype, weight_type( m_basis_grid, utils::ublas::project( m_basis_grid , weif::massdimm_weight( mass_response.first, mass_response.second, dimm_response.first, dimm_response.second, inner_diam, outer_diam, dimmsize, dimmbase ) ) ) ) );
+	std::pair< storage_type::iterator, bool > ret = m_storage.insert( std::make_pair( spectype, weight_type( m_basis_grid, utils::ublas::project( m_basis_grid , atmos::weight::massdimm_weight( mass_response.first, mass_response.second, dimm_response.first, dimm_response.second, inner_diam, outer_diam, dimmsize, dimmbase ) ) ) ) );
 
 	return ret.first->second;
 }
+
+}}

atmos-2.96.9.tar.gz/src/spectral_factory.h -> atmos-2.97.1.tar.gz/src/spectral_factory.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: spectral_factory.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: spectral_factory.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -23,29 +23,31 @@
 #include <string>
 
 #include <boost/filesystem/path.hpp>
-#include <boost/numeric/ublas/vector.hpp>
 #include <boost/function.hpp>
+#include <boost/numeric/ublas/vector.hpp>
 
+#include "spline.h"
 #include "nocase_less.h"
 #include "vector_generator.h"
+#include "utils.h"
 #include "weif.h"
 
-namespace spectral {
+namespace atmos {
+namespace weight {
 
 using namespace boost::numeric::ublas;
-typedef vector< double > vector_type;
+using namespace utils::algo;
+using namespace utils::ublas;
 
 std::pair< vector_type, vector_type > load_reaction( const boost::filesystem::path& file );
 std::pair< vector_type, vector_type > photon_distribution( const vector_type& spectral_lambda, const vector_type& spectral_response, const vector_type& device_lambda, const vector_type& device_response );
 
-void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid = utils::ublas::vector_generator< lg_generator >( lg_generator(0.0, 0.04, 0.25), 50 ) );
-
-}
+void dump_weight( const boost::filesystem::path& file, const boost::function1< vector_type, double >& function, const vector_type& grid = vector_generator< lg_generator >( lg_generator(0.0, 0.04, 0.25), 50 ) );
 
-class spectral_factory{
+class factory{
 public:
-	typedef boost::numeric::ublas::vector< double > vector_type;
-	typedef boost::numeric::ublas::vector< boost::numeric::ublas::vector<double> > vector_vector_type;
+	typedef atmos::weight::vector_type vector_type;
+	typedef vector< vector_type > vector_vector_type;
 
 	typedef utils::ublas::cubic_spline< vector_type, vector_vector_type, vector_type, vector_vector_type > weight_type;
 	typedef std::map< std::string, weight_type, utils::algo::nocase_less< std::string > > storage_type;
@@ -56,7 +58,7 @@
 	boost::filesystem::path m_ccd_response;
 	boost::filesystem::path m_wdir;
 
-	utils::ublas::vector_generator< lg_generator > m_basis_grid;
+	vector_generator< lg_generator > m_basis_grid;
 
 	storage_type m_storage;
 
@@ -65,14 +67,15 @@
 private:
 	const boost::filesystem::path get_spectrum_filename( const std::string& spectype ) const;
 public:
-	spectral_factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response );
-	~spectral_factory();
+	factory(const boost::filesystem::path& wd, const boost::filesystem::path& mass_response, const boost::filesystem::path& ccd_response );
+	~factory();
 
-	const weight_type& spectral( const std::string& spectype, const boost::numeric::ublas::vector<double>& inner_diam, const boost::numeric::ublas::vector<double>& outer_diam, double dimmbase, double dimmsize );
+	const weight_type& weight( const std::string& spectype, const vector_type& inner_diam, const vector_type& outer_diam, double dimmbase, double dimmsize );
 
 	const const_iterator begin() const { return m_storage.begin(); }
 	const const_iterator end()   const { return m_storage.end(); }
 };
 
-#endif //_SPECTRAL_FACTORY_H
+}}
 
+#endif //_SPECTRAL_FACTORY_H

atmos-2.96.9.tar.gz/src/table_storage.h -> atmos-2.97.1.tar.gz/src/table_storage.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: table_storage.h 143 2010-04-04 08:42:34Z matwey $
+    $Id: table_storage.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -75,29 +75,97 @@
 		m_storage.erase( unstable_remove_if( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
 	}
 	template<class Predicate> void remove_if( Predicate pred ) {
-		m_storage.erase( remove( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
+		using std::remove_if;
+		m_storage.erase( remove_if( m_storage.begin(), m_storage.end(), pred ), m_storage.end() );
 	}
 };
 
-template<class T, class F> class table_storage_column {
+namespace detail {
+template<class Result, class Iterator> struct table_storage_column_traits {
+	typedef Iterator subiterator;
+	typedef typename std::iterator_traits< subiterator >::reference arg_type;
+	typedef Result result_type;
+	typedef boost::function1< result_type, arg_type > function_type;
+	typedef boost::transform_iterator< function_type, subiterator, result_type > iterator;
+};
+}
+
+template<class Storage, class T> class table_storage_column_getter{
 public:
-	typedef T                subvalue_type;
-	typedef table_storage<T> table_storage_type;
-	typedef F                proxy_function_type;
+	typedef Storage storage_type;
+	typedef typename storage_type::const_iterator subiterator;
+	typedef const T& result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+private:
+	const storage_type& m_storage;
+	function_type m_function;
+public:
+	table_storage_column_getter( const storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
 
-	typedef typename table_storage_type::const_iterator const_subiterator;
+	iterator begin() const { return iterator( m_storage.begin(), m_function ); }
+	iterator end() const   { return iterator( m_storage.end(), m_function ); }
+};
 
-	typedef boost::transform_iterator< proxy_function_type, const_subiterator > const_iterator;
+template<class Storage, class T> class table_storage_column_setter{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::iterator subiterator;
+	typedef T& result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
 private:
-	const table_storage<T>& m_storage;
-	proxy_function_type m_proxy_function;
+	storage_type& m_storage;
+	function_type m_function;
 public:
-	table_storage_column( const table_storage<T>& storage, const F& proxy_function ):
-		m_storage( storage ), m_proxy_function( proxy_function ) {
-	}
+	table_storage_column_setter( storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
+
+	iterator begin() { return iterator( m_storage.begin(), m_function ); }
+	iterator end()   { return iterator( m_storage.end(), m_function ); }
+};
 
-	const_iterator begin() const { return const_iterator( m_storage.begin(), m_proxy_function ); }
-	const_iterator end()   const { return const_iterator( m_storage.end(), m_proxy_function ); }
+template<class Storage, class T> class table_storage_column_generator{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::const_iterator subiterator;
+	typedef T result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+private:
+	const storage_type& m_storage;
+	function_type m_function;
+public:
+	table_storage_column_generator( const storage_type& storage, const function_type& function ): m_storage(storage), m_function(function) {}
+
+	iterator begin() const { return iterator( m_storage.begin(), m_function ); }
+	iterator end() const   { return iterator( m_storage.end(), m_function ); }
+};
+
+template<class Storage, class T> class table_storage_column_proxy{
+public:
+	typedef Storage storage_type;
+	typedef typename storage_type::iterator       subiterator;
+	typedef typename storage_type::const_iterator const_subiterator;
+	typedef T& result_type;
+	typedef const T& const_result_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::function_type function_type;
+	typedef typename detail::table_storage_column_traits< result_type, subiterator >::iterator      iterator;
+	typedef typename detail::table_storage_column_traits< const_result_type, const_subiterator >::function_type const_function_type;
+	typedef typename detail::table_storage_column_traits< const_result_type, const_subiterator >::iterator      const_iterator;
+private:
+	storage_type& m_storage;
+	function_type m_function;
+	const_function_type m_const_function;
+public:
+	table_storage_column_proxy( storage_type& storage, const function_type& function, const const_function_type& const_function ):
+		m_storage(storage), m_function(function), m_const_function(const_function) {}
+	template<class F> table_storage_column_proxy( storage_type& storage, const F& function ):
+		m_storage(storage), m_function(function), m_const_function(function) {}
+
+	iterator begin() { return iterator( m_storage.begin(), m_function ); }
+	iterator end()   { return iterator( m_storage.end(), m_function ); }
+	const_iterator begin() const { return const_iterator( m_storage.begin(), m_const_function ); }
+	const_iterator end()   const { return const_iterator( m_storage.end(), m_const_function ); }
 };
 
 template<class T> void swap( table_storage<T>& t, table_storage<T>& u ) {

atmos-2.96.9.tar.gz/src/typesdef.h -> atmos-2.97.1.tar.gz/src/typesdef.h Changed

atmos-2.97.1.tar.gz/src/utils.cpp Added

@@ -0,0 +1,48 @@
+/*
+    $Id: utils.cpp 180 2011-01-30 15:56:33Z matwey $
+    Copyright (C) 2011  Sternberg Astronomical Institute, MSU
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <cmath>
+
+#include "utils.h"
+
+namespace utils {
+namespace algo {
+
+linear_generator::linear_generator( result_type start, result_type step ): start_(start), step_(step) {
+}
+linear_generator::result_type linear_generator::operator()( size_type i ) const {
+	return ( start_ + step_ * i );
+}
+linear_generator::size_type linear_generator::size( result_type max ) const {
+	return static_cast<size_type>( std::floor( (max - start_) / step_ ) );
+}
+
+lg_generator::lg_generator( result_type start, result_type lstep, result_type gstep ): start_(start), t_(lstep), q_(result_type(1) + gstep / 2) {
+	max_ = static_cast< size_type >( std::floor( result_type(2) / gstep - start / lstep ) );
+}
+lg_generator::result_type lg_generator::operator()( size_type i ) const {
+	if( i <= max_ ) return ( start_ + t_ * i );
+	return ( start_ + t_ * max_ ) * std::pow( q_, static_cast<int>(i - max_) ); // FIXME: pow(double,int) for <4.2.1
+}
+lg_generator::size_type lg_generator::size( result_type max ) const {
+	result_type gmax = std::floor( std::log( max / ( start_ + t_ * max_ ) ) / std::log( q_ ) );
+	if( gmax > result_type(0) ) return static_cast< size_type >( gmax ) + max_ + 1;
+	return static_cast< size_type >(std::floor( (max - start_) / t_ ));
+}
+
+}}

atmos-2.96.9.tar.gz/src/utils.h -> atmos-2.97.1.tar.gz/src/utils.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: utils.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: utils.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -16,6 +16,13 @@
     along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
+#ifndef _UTILS_H_
+#define _UTILS_H_
+
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/triangular.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
 namespace utils {
 namespace ublas {
 
@@ -60,5 +67,31 @@
 	return ret;
 }
 
-} // utils
 } // ublas
+
+namespace algo {
+struct linear_generator {
+	typedef double result_type;
+	typedef size_t size_type;
+
+	linear_generator( result_type start, result_type step );
+	result_type operator()( size_type i ) const;
+	size_type size( result_type max ) const;
+private:
+	result_type start_,step_;
+};
+struct lg_generator {
+	typedef double result_type;
+	typedef size_t size_type;
+
+	lg_generator( result_type start, result_type lstep, result_type gstep );
+	result_type operator()( size_type i ) const;
+	size_type size( result_type max ) const;
+private:
+	result_type start_,t_,q_;
+	size_type max_;
+};
+} // algo
+} // utils
+
+#endif // _UTILS_H_

atmos-2.96.9.tar.gz/src/weif.cpp -> atmos-2.97.1.tar.gz/src/weif.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: weif.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: weif.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -17,16 +17,18 @@
 */
 
 #include <cmath>
-#include <boost/function.hpp>
-
-#include <boost/numeric/ublas/io.hpp>
 
-#include "weif.h"
+#include <boost/function.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
 
-#include "gsl_wrap.h"
+#include "atmos.h"
 #include "fft.h"
+#include "gsl_wrap.h"
+#include "utils.h"
+#include "weif.h"
 
-namespace weif {
+namespace atmos {
+namespace weight {
 
 const double POWIND = -2.66666666667;
 const double POWDIFF = -0.666666666667;
@@ -45,8 +47,9 @@
 	double diameter3;
 	double diameter4;
 	double PIx2lambda0;
-	const sfun_type& sfun;
+	sfun_type sfun;
 	double altitude;
+	double powind;
 
 	mass_params( const sfun_type& sfun ): sfun(sfun) {}
 };
@@ -58,31 +61,29 @@
 	double dimmbase;
 	double PIx2lambda0;
 	double phi;
-	const sfun_type& sfun;
+	sfun_type sfun;
 	double altitude;
 
 	dimm_params( const sfun_type& sfun ): sfun(sfun) {}
 };
 
 double mass_integrand(double x, void * params) {
+	using namespace std;
+	
 	mass_params p = *( mass_params *) params;
-	double apert1 = 0.0, apert2 = 0.0, apert3 = 0.0, apert4 = 0;
+	double p1 = 0;
+	double p3 = 0;
 	double f = p.altitude*x*x/2;
 	const std::complex<double>& sfun_f = p.sfun(f);
 	double sf = imag( sfun_f ) * cos( p.PIx2lambda0*f ) + real( sfun_f ) * sin( p.PIx2lambda0*f );
-	sf = sf*sf*pow( x, POWIND );
-	if( p.diameter1 != p.diameter2 ) {                // First aperture:
-		apert1 = p.diameter1*j1( p.diameter1*x );
-		apert2 = p.diameter2*j1( p.diameter2*x );
-		apert1 = ( apert1 - apert2 )/( p.diameter1*p.diameter1 - p.diameter2*p.diameter2 );
+	sf = sf*sf*pow( x, p.powind );
+	if( p.diameter1 != p.diameter2 ) {
+		p1 = ( p.diameter1*j1( p.diameter1*x ) - p.diameter2*j1( p.diameter2*x ) ) / ( p.diameter1*p.diameter1 - p.diameter2*p.diameter2 );
 	}
-	if( p.diameter3 != p.diameter4) {                 // Second aperture: the same with indices 1->3, 2->4:
-		apert3 = p.diameter3*j1( p.diameter3*x );
-		apert4 = ( p.diameter4 == 0 )? 0 : p.diameter4*j1( p.diameter4*x );
-		apert3 = ( apert3 - apert4 )/( p.diameter3*p.diameter3 - p.diameter4*p.diameter4 );
+	if( p.diameter3 != p.diameter4 ) {
+		p3 = ( p.diameter3*j1( p.diameter3*x ) - p.diameter4*j1( p.diameter4*x ) ) / ( p.diameter3*p.diameter3 - p.diameter4*p.diameter4 );
 	}
-	apert1 = 2.0*( apert1 - apert3 )/x;
-	return sf*apert1*apert1;
+	return sf * 4.0 * (p1 * p3) / (x*x);
 }
 
 double dimm_integrand(double x, void * params) {
@@ -113,12 +114,12 @@
 	return it;
 }
 
-std::pair< vector<double>, vector< std::complex<double> > > spectral_band_fourier( const vector<double>& lambda, const vector<double>& dist, unsigned int center ) {
+std::pair< vector_type, vector< complex_type > > spectral_band_fourier( const vector<double>& lambda, const vector<double>& dist, unsigned int center ) {
 	const int GRIDFACT = 4 ;    // 0 - no stretch, 1 - stretch twice, 2 - stretch 4-times etc.
 	double dlambda = lambda(1) - lambda(0);
 	int dist_ftnn = 2*GRIDFACT*nexthigher(dist.size());
 	vector< std::complex<double> > dist_ft = zero_vector< std::complex<double> >( dist_ftnn );
-	vector<double> dist_freq( dist_ftnn/2 );
+	vector< double > dist_freq( dist_ftnn/2 );
 	for( vector<double>::size_type l = 0; l < dist_freq.size(); l++ ) dist_freq(l) = l/(dlambda*MKM2CM)/dist_ftnn;
 	for( vector<double>::size_type l = 0; l < dist.size(); l++ ) {      // Fill in the data array for FFT:
 		vector<double>::size_type j = ( l < center ) ? dist_ftnn + l - center : l - center;
@@ -128,10 +129,12 @@
 	return std::make_pair( dist_freq, subrange( dist_ft, 0, dist_ftnn/2 ) );
 }
 
-double mass_weight( double height, double eps1, double eps2, double eps3, double eps4, double lambda0, const boost::function1< std::complex<double>, double >& sfun, double max_freq ) {
+namespace detail {
+
+basic_weight::value_type basic_weight::do_calculate( value_type height, value_type eps1, value_type eps2, value_type eps3, value_type eps4 ) const {
 	gsl::integration_workspace w(WORKSPACE);
 	double error;
-	mass_params params(sfun);
+	mass_params params( mass_sfourier_ );
 	gsl_function F;
 	F.function = &mass_integrand;
 	F.params = &params;
@@ -139,31 +142,60 @@
 	params.diameter2 = M_PI*eps2;
 	params.diameter3 = M_PI*eps3;
 	params.diameter4 = M_PI*eps4;
-	params.PIx2lambda0 = 2*M_PI*lambda0*MKM2CM;
+	params.PIx2lambda0 = 2*M_PI * mass_leff_ * MKM2CM;
+	params.powind = POWIND + freq_power_;
 	double result;
 	params.altitude = height*KM2CM;
-	double hi_limit = std::sqrt( 2 * max_freq / params.altitude );
+	double hi_limit = std::sqrt( 2 * mass_freq_max_ / params.altitude );
 	if( std::isfinite( hi_limit ) )
 		gsl::integration_qags( &F, 0, hi_limit, 0, EPS, WORKSPACE, w, &result, &error);
 	else
 		gsl::integration_qagiu( &F, 0, 0, EPS, WORKSPACE, w, &result, &error);
-	return result*CALIBR;
+	return result*CALIBR*std::pow(M_PI*100,freq_power_);
+}
+
+basic_weight::basic_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& freq_power ):
+	ap_inn_( ap_inn ),
+	ap_out_( ap_out ),
+	mass_leff_index_( lambda_effective( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ) ).index() ),
+	mass_leff_( mass_lambda( mass_leff_index_ ) ),
+	mass_sfourier_( spectral_band_fourier( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ), mass_leff_index_ ) ),
+	freq_power_(freq_power) {
+
+	mass_freq_max_ = mass_sfourier_.first_data()( mass_sfourier_.size() - 1 );
+}
+
+basic_weight::result_type basic_weight::operator() ( const arg_type& height ) const {
+	const size_type mass_nchan = atmos::channels;
+	vector_type ret( atmos::indices_size );
+	size_type cur_n = 0;
+
+	for( size_type i = 0; i < mass_nchan; ++i )
+		ret( cur_n++ ) = do_calculate( height, ap_out_(i), ap_inn_(i), ap_out_(i), ap_inn_(i) );
+
+	for( size_type i = 0; i < mass_nchan; ++i )
+		for( size_type j = i + 1; j < mass_nchan; j++ )
+			ret( cur_n++ ) = do_calculate( height, ap_out_(i), ap_inn_(i), ap_out_(j), ap_inn_(j) );
+
+	return ret;
+}
+
 }
 
-double dimm_weight( double height, double diam, double dbase, double phi, double lambda0, const boost::function1< std::complex<double>, double >& sfun, double max_freq ) {
+massdimm_weight::value_type massdimm_weight::do_dimm_calculate( value_type height, value_type diam, value_type dbase, value_type phi ) const {
 	gsl::integration_workspace w(WORKSPACE);
 	double error;
-	dimm_params params(sfun);
+	dimm_params params(dimm_sfourier_);
 	gsl_function F;
 	F.function = &dimm_integrand;
 	F.params = &params;
 	params.aperture = M_PI*diam;
 	params.dimmbase = 2*M_PI*dbase;
-	params.PIx2lambda0 = 2*M_PI*lambda0*MKM2CM;
+	params.PIx2lambda0 = 2*M_PI*dimm_leff_*MKM2CM;
 	params.phi = phi;
 	double result;
 	params.altitude = height*KM2CM;
-	double hi_limit = std::sqrt( 2 * max_freq / params.altitude );
+	double hi_limit = std::sqrt( 2 * dimm_freq_max_ / params.altitude );
 	if( std::isfinite( hi_limit ) )
 		gsl::integration_qags( &F, 0, hi_limit, 0, EPS, WORKSPACE, w, &result, &error);
 	else
@@ -171,57 +203,26 @@
 	return result*DIMMCONST;
 }
 
-}
-
-weif::massdimm_weight::massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, double dimmsize, double dimmbase ):
-	ap_inn_( ap_inn ), ap_out_( ap_out ), dimmsize_( dimmsize * 100 ), dimmbase_( dimmbase * 100 ),
-	mass_leff_index_( lambda_effective( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ) ).index() ),
+massdimm_weight::massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& dimmsize, const value_type& dimmbase ):
+	basic_weight( mass_lambda, mass_response, ap_inn, ap_out ),
+	dimmsize_( dimmsize * 100 ), dimmbase_( dimmbase * 100 ),
 	dimm_leff_index_( lambda_effective( dimm_lambda, dimm_response ).index() ),
-	mass_leff_( mass_lambda( mass_leff_index_ ) ),
 	dimm_leff_( dimm_lambda( dimm_leff_index_ ) ),
-	mass_sfourier_( spectral_band_fourier( mass_lambda, element_div( mass_response, mass_lambda*MKM2CM ), mass_leff_index_ ) ),

atmos-2.96.9.tar.gz/src/weif.h -> atmos-2.97.1.tar.gz/src/weif.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: weif.h 158 2010-09-19 08:38:25Z matwey $
+    $Id: weif.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2010  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,70 +19,66 @@
 #ifndef _WEIF_H
 #define _WEIF_H
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/vector_proxy.hpp>
+#include <complex>
 
-#include <utility>
+#include <boost/numeric/ublas/vector.hpp>
 
 #include "spline.h"
 
-namespace weif {
-
-	using namespace std;
+namespace atmos {
+namespace weight {
 	using namespace boost::numeric::ublas;
-	using boost::numeric::ublas::vector;
 
-	class massdimm_weight {
+	typedef double value_type;
+	typedef std::complex< value_type > complex_type;
+	typedef vector< value_type >   vector_type;
+	typedef vector_type::size_type size_type;
+	typedef vector_type            weight_type;
+	
+namespace detail {
+	class basic_weight {
 	public:
-		typedef double value_type;
-		typedef vector< value_type > vector_type;
-		typedef vector_type::size_type size_type;
-		typedef vector_type result_type;
-
-		typedef utils::ublas::cubic_spline< vector_type, vector< std::complex<double> >, vector_type, vector< std::complex<double> > > sfourier_type;
+		typedef atmos::weight::value_type value_type;
+		typedef value_type  arg_type;
+		typedef weight_type result_type;
+
+	protected:
+		typedef utils::ublas::cubic_spline< vector_type, vector< complex_type >, vector_type, vector< complex_type > > sfourier_spline_type;
+		typedef sfourier_spline_type sfourier_type;
+	private:
+		value_type do_calculate( value_type height, value_type eps1, value_type eps2, value_type eps3, value_type eps4 ) const;
 	public:
-		massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, double dimmsize, double dimmbase );
+		basic_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& freq_power = value_type(0) );
 
-		vector_type operator() ( const value_type& height ) const;
+		virtual result_type operator() ( const arg_type& height ) const;
 	private:
 		vector_type ap_inn_;
 		vector_type ap_out_;
-		double dimmsize_;
-		double dimmbase_;
 		size_type mass_leff_index_;
-		size_type dimm_leff_index_;
-		double mass_leff_;
-		double dimm_leff_;
+		value_type mass_leff_;
 		sfourier_type mass_sfourier_;
-		sfourier_type dimm_sfourier_;
-		double mass_freq_max_;
-		double dimm_freq_max_;
+		value_type mass_freq_max_;
+		value_type freq_power_;
 	};
-
-	/* make_weight( crv, geometry ) */
 }
 
-    struct linear_generator {
-        typedef double result_type;
-        typedef size_t size_type;
-
-        linear_generator( result_type start, result_type step );
-        result_type operator()( size_type i ) const;
-        size_type size( result_type max ) const;
-    private:
-        result_type start_,step_;
-    };
-
-    struct lg_generator {
-        typedef double result_type;
-        typedef size_t size_type;
-
-        lg_generator( result_type start, result_type lstep, result_type gstep );
-        result_type operator()( size_type i ) const;
-        size_type size( result_type max ) const;
-    private:
-        result_type start_,t_,q_;
-        size_type max_;
-    };
+class massdimm_weight:
+	public detail::basic_weight {
+private:
+	value_type do_dimm_calculate( value_type height, value_type diam, value_type dbase, value_type phi ) const;
+public:
+	massdimm_weight( const vector_type& mass_lambda, const vector_type& mass_response, const vector_type& dimm_lambda, const vector_type& dimm_response, const vector_type& ap_inn, const vector_type& ap_out, const value_type& dimmsize, const value_type& dimmbase );
+	
+	result_type operator() ( const arg_type& height ) const;
+private:
+	value_type dimmsize_;
+	value_type dimmbase_;
+	size_type dimm_leff_index_;
+	value_type dimm_leff_;
+	sfourier_type dimm_sfourier_;
+	value_type dimm_freq_max_;
+};
+
+}}
 
 #endif // _WEIF_H

atmos-2.96.9.tar.gz/src/worksheet.cpp -> atmos-2.97.1.tar.gz/src/worksheet.cpp Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: worksheet.cpp 158 2010-09-19 08:38:25Z matwey $
+    $Id: worksheet.cpp 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -22,11 +22,9 @@
 #include <boost/algorithm/string.hpp>
 #include <boost/bind.hpp>
 #include <boost/filesystem/path.hpp>
-#include <boost/iterator.hpp>
-#include <boost/iterator/filter_iterator.hpp>
 #include <boost/format.hpp>
 #include <boost/lexical_cast.hpp>
-#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
 
 #include <yalpa/sun.h>
 #include <yalpa/sidtime.h>
@@ -38,7 +36,7 @@
 #include "weif.h" // for wf::*_generator
 #include "worksheet.h"
 
-#include "photometry.h"
+namespace atmos {
 
 struct worksheet::profile_grid {
 	boost::numeric::ublas::vector<double> vec;
@@ -50,7 +48,7 @@
 	public local_context< worksheet > {
 public:
 	typedef context_type::input::time_type time_type;
-	typedef context_type::spectral_factory_type::weight_type weight_type;
+	typedef context_type::weight_factory_type::weight_type weight_type;
 private:
 	std::string m_spectral;
 	yalpa::equatorial<double> m_coords;
@@ -58,7 +56,7 @@
 	object( const yalpa::equatorial<double> ecoords, const std::string& spectral, context_reference_type ctx ): local_context_type( ctx ), m_spectral(spectral), m_coords(ecoords) {}
 
 	double airmass( const time_type& time ) const;
-	const weight_type& spectral();
+	const weight_type& weight();
 	const std::string& sclass() const { return m_spectral; }
 };
 
@@ -68,27 +66,29 @@
 class worksheet::input::background_accumulator:
 	public measurement_accumulator {
 private:
-	std::auto_ptr<background_model> m_model;
-	const ublas::vector< double > m_deadtime;
-	background_nonpoisson_filter m_nonpoisson_filter;
+	std::auto_ptr< atmos::background::model::abstract > m_model;
+	const atmos::parameters::vector_type m_deadtime;
+	atmos::background::filter::nonpoisson m_nonpoisson_filter;
 public:
-	background_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ):
+	background_accumulator( const time_type& time, context_reference_type ctx ):
 		measurement_accumulator( time, ctx ),
-		m_model( ctx.m_bfactory.create_new_model(time) ),
-		m_deadtime( ctx.m_pproxy.nonlinearity() ),
+		m_model( ctx.m_background_factory.create_new_model(time) ),
+		m_deadtime( ctx.parameters().nonlinearity() ),
 		m_nonpoisson_filter(
 			ctx.m_va"background-threshold".as<double>(),
-			ctx.pproxy().nonpoisson()(3),
-			1000 * ctx.pproxy().basetime() / ctx.pproxy().exposition() ) {
+			ctx.parameters().nonpoisson()(3),
+			1000 * ctx.parameters().basetime() / ctx.parameters().exposition() ) {
 
 		if( !m_model.get() )
 			throw std::logic_error("There is not any background model");
 	};
 	virtual ~background_accumulator() {};
 	virtual void accumulate(  const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+		using atmos::indices::flux_nonlinearity_correction;
+		
 		if( ! m_nonpoisson_filter( time, mean, vars ) ) {
 			if( m_model.get() ) {
-				m_model->accumulate( time, indices::flux_nonlinearity_correction( mean, m_deadtime ) );
+				m_model->accumulate( time, flux_nonlinearity_correction( mean, m_deadtime ) );
 			}
 		}
 	}
@@ -98,8 +98,8 @@
 		context().m_file_output.background( starttime() );
 
 		if( m_model.get() && m_model->size() > 0 ) {
-			context().m_bfactory.handle_new_model( m_model.release() );
-			context().m_file_output.fluxes( starttime(), context().m_bfactory.background( starttime() ) );
+			context().m_background_factory.handle_new_model( m_model.release() );
+			context().m_file_output.fluxes( starttime(), context().m_background_factory.background( starttime() ) );
 		}
 	}
 };
@@ -108,33 +108,35 @@
 public:
 	typedef context_type::input::size_type size_type;
 private:
-	boost::posix_time::ptime m_accumulation_begin;
-	boost::posix_time::ptime m_accumulation_end;
+	time_type m_accumulation_begin;
+	time_type m_accumulation_end;
 	double m_mean_mz;
 	size_type m_num;
 
-	dimm_storage m_dimm;
-
+	atmos::dimm::storage    m_dimm;
+	atmos::indices::storage m_indices;
+	
 	//FIXME: use params proxy lazy initialization
-	const ublas::vector< double > m_deadtime;
-	const ublas::vector< double > m_scatter;
-	const ublas::vector< double > m_npoisson;
+	const atmos::parameters::vector_type m_deadtime;
+	const atmos::parameters::vector_type m_scatter;
+	const atmos::parameters::vector_type m_npoisson;
 
-	indices_storage m_indices;
-	indices_flux_threshold_filter m_flux_filter;
+	atmos::indices::filter::flux_threshold m_flux_filter;
 public:
-	normal_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ):
+	normal_accumulator( const time_type& time, context_reference_type ctx ):
 		measurement_accumulator( time, ctx ),
 		m_mean_mz(0.0),
 		m_num(0),
-		m_deadtime( ctx.pproxy().nonlinearity() ),
-		m_scatter( ctx.pproxy().scatter() ),
-		m_npoisson( ctx.pproxy().nonpoisson() ),
+		m_deadtime( ctx.parameters().nonlinearity() ),
+		m_scatter( ctx.parameters().scatter() ),
+		m_npoisson( ctx.parameters().nonpoisson() ),
 		m_flux_filter( ctx.m_va"flux-threshold".as<double>() ) {
 
 	}
 	virtual ~normal_accumulator() {};
 	virtual void accumulate( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {
+		using namespace atmos::indices;
+		
 		if( m_accumulation_begin == boost::posix_time::not_a_date_time )
 			m_accumulation_begin = time;
 		m_accumulation_end = time;
@@ -142,38 +144,42 @@
 		m_num++;
 		m_mean_mz = m_mean_mz * ( double(m_num-1) / double(m_num) ) + context().curobject().airmass( time ) * ( 1.0/m_num );
 
-		indices::flux_type flux_corrected = indices::flux_nonlinearity_correction( mean, m_deadtime );
+		flux_type flux_corrected = flux_nonlinearity_correction( mean, m_deadtime );
 
-		const boost::tuple< vars_type, rho1_type, rho2_type >& vars_corrected = indices::vars_correction(
+		const boost::tuple< vars_type, rho1_type, rho2_type >& vars_corrected = vars_correction(
 			flux_corrected,
 			element_prod( vars, outer_prod( mean, mean ) ),
 			element_prod( rho1, outer_prod( mean, mean ) ),
 			element_prod( rho2, element_prod( mean, mean ) ), m_deadtime, m_npoisson );
 
-		flux_corrected = indices::flux_background_correction( flux_corrected, m_scatter, context().m_bfactory.background( time ) );
+		flux_corrected = flux_background_correction( flux_corrected, m_scatter, context().background_factory().background( time ) );
 
 		if( ! m_flux_filter( flux_corrected ) ) {
-			const boost::tuple< indices::indices_type, indices::indesis_type>& inds = indices::calculate( flux_corrected, vars_corrected.get<0>(), vars_corrected.get<1>(), vars_corrected.get<2>() );
+			const boost::tuple< indices_type, indesis_type>& inds = calculate( flux_corrected, vars_corrected.get<0>(), vars_corrected.get<1>(), vars_corrected.get<2>() );
 
 			m_indices.push_back( time, flux_corrected, inds.get<0>(), inds.get<1>() );
 		}
 	}
 	virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) {
+		using atmos::dimm::correction;
+		
 		if( m_accumulation_begin == boost::posix_time::not_a_date_time )
 			m_accumulation_begin = time;
 		m_accumulation_end = time;
 
-		m_dimm.push_back( time, framenum, left, right, dimm::correction( xvars ), dimm::correction( yvars ) );
+		m_dimm.push_back( time, framenum, left, right, correction( xvars ), correction( yvars ) );
 	}
 	virtual void flush() {
 		using namespace utils::algo;
+		using namespace atmos;
+		using namespace boost::numeric::ublas;
 
 		const boost::posix_time::ptime& meantime = m_accumulation_begin + (m_accumulation_end-m_accumulation_begin)/2;
 
 		typedef context_type::integrals_factory_type integrals_factory_type;
 		typedef integrals_factory_type::functor_type integrals_functor_type;
 		typedef integrals_functor_type::result_type  integrals_result_type;
-		typedef restoration_functor                   restoration_functor_type;
+		typedef profile::functor                    restoration_functor_type;
 		typedef restoration_functor_type::result_type restoration_result_type;
 
 		context().m_file_output.normal( starttime(), m_mean_mz );
@@ -185,12 +191,17 @@
 			double dimm_turbulence_err = 0.0;
 
 			if( m_dimm.size() > 1 ) {
-				m_dimm.correct_low_vars();

atmos-2.96.9.tar.gz/src/worksheet.h -> atmos-2.97.1.tar.gz/src/worksheet.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: worksheet.h 147 2010-04-17 09:25:25Z matwey $
+    $Id: worksheet.h 180 2011-01-30 15:56:33Z matwey $
     Copyright (C) 2009  Sternberg Astronomical Institute, MSU
 
     This program is free software: you can redistribute it and/or modify
@@ -19,58 +19,45 @@
 #ifndef _WORKSHEET_H
 #define _WORKSHEET_H
 
+#include <boost/program_options.hpp>
+
 #include "background_factory.h" // for background_factory
 #include "file_input.h"       // for file_input
 #include "file_output.h"      // for file_output
-#include "params.h"           // for parameters_*
-#include "spectral_factory.h" // for spectral_factory
-#include "typesdef.h"         // for nocase_less
-
 #include "indices.h"
 #include "integrals.h"
+#include "local_context.h"
+#include "params.h"           // for parameters_*
 #include "profile.h"
+#include "spectral_factory.h" // for spectral_factory
 
-#include <boost/date_time/posix_time/posix_time.hpp>
-#include <boost/program_options.hpp>
+namespace atmos {
 
-template<class Context> class local_context {
-	public:
-		typedef Context                context_type;
-		typedef context_type&          context_reference_type;
-		typedef const context_type&    context_const_reference_type;
-		typedef local_context<Context> local_context_type;
-	private:
-		context_reference_type m_ctx;
-	protected:
-		inline context_reference_type       context()       throw() { return m_ctx; }
-		inline context_const_reference_type context() const throw() { return m_ctx; }
-	public:
-		local_context( context_reference_type ctx ) throw() : m_ctx(ctx) {}
-};
+using utils::local_context;
 
 class worksheet {
 	friend class local_context<worksheet>;
 	public:
 		typedef local_context<worksheet> context_type;
-		typedef parameters_container< std::map< std::string, std::string, utils::algo::nocase_less< std::string > > > parameters_container_type;
-		typedef parameters_proxy< parameters_container_type > parameters_proxy_type;
-		typedef spectral_factory spectral_factory_type;
-		typedef background_factory background_factory_type;
-		typedef integrals_factory integrals_factory_type;
+		typedef atmos::parameters::container      parameters_container_type;
+		typedef atmos::parameters::proxy          parameters_proxy_type;
+		typedef atmos::weight::factory            weight_factory_type;
+		typedef atmos::background::factory        background_factory_type;
+		typedef atmos::integrals::functor_factory integrals_factory_type;
 	public:
 		static int parse_program_options( int argc, char** argv, boost::program_options::variables_map& va );
 	public:
 		struct profile_grid;
 	private:
 		class input:
-			public file_input,
+			public atmos::input::file,
 			public local_context< worksheet > {
 		private:
 			class measurement_accumulator: public local_context< worksheet > {
 			private:
-				boost::posix_time::ptime m_starttime;
+				time_type m_starttime;
 			public:
-				measurement_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ) throw() : local_context_type( ctx ), m_starttime(time) {};
+				measurement_accumulator( const time_type& time, context_reference_type ctx ) throw() : local_context_type( ctx ), m_starttime(time) {};
 				virtual ~measurement_accumulator() = 0;
 				virtual void accumulate( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) = 0;
 				virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) = 0;
@@ -80,12 +67,12 @@
 					return dynamic_cast<const T*>(this);
 				}
 
-				inline const boost::posix_time::ptime& starttime() const { return m_starttime; }
+				inline const time_type& starttime() const { return m_starttime; }
 			};
 
 			class dummy_accumulator: public measurement_accumulator {
 			public:
-				dummy_accumulator( const boost::posix_time::ptime& time, context_reference_type ctx ) throw(): measurement_accumulator( time, ctx ) {};
+				dummy_accumulator( const time_type& time, context_reference_type ctx ) throw(): measurement_accumulator( time, ctx ) {};
 				virtual ~dummy_accumulator() {};
 				virtual void accumulate(  const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 ) {}
 				virtual void dimm_accumulate( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars ) {}
@@ -97,7 +84,7 @@
 		private:
 			std::auto_ptr<measurement_accumulator> m_maccum;
 
-			template<class T> void switch_accumulator( const boost::posix_time::ptime& time ) {
+			template<class T> void switch_accumulator( const time_type& time ) {
 				std::auto_ptr< measurement_accumulator > new_accumulator( new T( time,context() ) );
 				std::swap( m_maccum, new_accumulator );
 				new_accumulator->flush();
@@ -116,7 +103,7 @@
 			void on_measurements( const time_type& time, const flux_type& mean, const vars_type& vars, const rho1_type& rho1, const rho2_type& rho2 );
 			void on_dimm( const time_type& time, unsigned int framenum, const dimm_flux_type& left, const dimm_flux_type& right, const dimm_vars_type& xvars, const dimm_flux_type& yvars );
 		public:
-			input( std::istream& stm, context_reference_type ctx ): file_input( stm ), local_context_type( ctx ), m_maccum( new dummy_accumulator( time_type() ,ctx) ) {}
+			input( std::istream& stm, context_reference_type ctx ): atmos::input::file( stm ), local_context_type( ctx ), m_maccum( new dummy_accumulator( time_type() ,ctx) ) {}
 		};
 		
 		class object;
@@ -125,29 +112,26 @@
 		
 		parameters_container_type m_pcontainer;
 		parameters_proxy_type     m_pproxy;
-
-		spectral_factory_type     m_sfactory;
-
-		background_factory_type   m_bfactory;
-
-		integrals_factory_type    m_ifactory;
+		weight_factory_type       m_weight_factory;
+		background_factory_type   m_background_factory;
+		integrals_factory_type    m_integrals_factory;
 		
-		input       m_file_input;
-		file_output m_file_output;
+		input        m_file_input;
+		output::file m_file_output;
 
 		std::auto_ptr<object> m_curobject;
 	private:
-		background_model* select_background_model( const boost::posix_time::ptime& time ) const;
+		background::model::abstract* select_background_model( const background::time_type& time ) const;
 	public:
 		worksheet( const boost::program_options::variables_map& va );
 		~worksheet( );
 
-		inline const parameters_proxy_type& pproxy()     const { return m_pproxy;   }
-		inline const spectral_factory_type& sfactory()   const { return m_sfactory; }
-		inline const background_factory_type& bfactory() const { return m_bfactory; }
-		inline const integrals_factory_type& ifactory()  const { return m_ifactory; }
+		inline const parameters_proxy_type&   parameters()         const { return m_pproxy;   }
+		inline const weight_factory_type&     weight_factory()     const { return m_weight_factory; }
+		inline const background_factory_type& background_factory() const { return m_background_factory; }
+		inline const integrals_factory_type&  integrals_factory()  const { return m_integrals_factory; }
 
-		double sun_altitude( const boost::posix_time::ptime& time ) const;
+		double sun_altitude( const background::time_type& time ) const;
 
 		inline object& curobject() {
 			if( ! m_curobject.get() )
@@ -160,4 +144,6 @@
 		void post();
 };
 
+}
+
 #endif // _WORKSHEET_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/Makefile.am -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/Makefile.am Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/Makefile.in -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/Makefile.in Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/bidiagonal_transform.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/bidiagonal_transform.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/cholesky_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/cholesky_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: cholesky_decomposition.h 68 2009-04-20 18:46:20Z matwey $
+    $Id: cholesky_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,24 +19,23 @@
 #ifndef _CHOLESKY_DECOMPOSITION_H
 #define _CHOLESKY_DECOMPOSITION_H
 
-#include <limits>
-
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/banded.hpp>
 #include <boost/numeric/ublas/triangular.hpp>
 #include <boost/numeric/ublas/vector_proxy.hpp>
 
-using namespace boost::numeric::ublas;
+#include <limits>
 
 namespace lsp{
-	
-	
+
+using namespace boost::numeric::ublas;
+
 /**
  *  @brief Upper Cholesky decomposition
  *  @paramin,out m The matrix to be decomposited. Result is stored in this object.
  *
- *  Cholesky decomposition is defined only for symmetric positive semidefined matrixes.
+ *  Cholesky decomposition is defined only for symmetric positive semidefinite matrixes.
  *  This conditions aren't checked.
  *  
  *  The decomposition defined as \f$ C = R^T R \f$ where \f$ R \f$ is upper triangular.

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/givens_rotation.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/givens_rotation.h Changed

@@ -1,6 +1,6 @@
 /*
-    $Id: givens_rotation.h 69 2009-08-28 18:32:24Z matwey $
-    Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+    $Id: givens_rotation.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2008,2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -19,14 +19,14 @@
 #ifndef _GIVENS_ROTATION_H
 #define _GIVENS_ROTATION_H
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 
-using namespace boost::numeric::ublas;
+#include <lsp/utils.h>
 
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class givens_rotation
  *  @brief A functor for the Givens rotation transformation
@@ -133,7 +133,7 @@
  *
  */
 	template<class U> void apply ( U& x, U& y ) const {
-		U w ( x * m_c + y * m_s );
+		typename lsp::temporary_type_traits< U >::type w ( x * m_c + y * m_s );
 		y = x * ( -m_s ) + y * m_c;
 		x = w;
 	}

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/householder_transform.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/householder_transform.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: householder_transform.h 60 2009-03-01 18:26:41Z matwey $
+    $Id: householder_transform.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,18 +19,19 @@
 #ifndef _HOUSEHOLDER_TRANSFORM_H
 #define _HOUSEHOLDER_TRANSFORM_H
 
-#include <lsp/utils.h>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/symmetric.hpp>
 
 #include <algorithm>
 #include <limits>
 
-#include <boost/numeric/ublas/vector.hpp>
-#include <boost/numeric/ublas/matrix.hpp>
-
-using namespace boost::numeric::ublas;
+#include <lsp/utils.h>
 
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class householder_transform
  *  @brief A functor for the Householder transformation
@@ -100,7 +101,7 @@
 			m_s += ( v(p)/w )*( v(p)/w );
 			for( size_type i = l; i < m; ++i )
 				m_s += ( v(i)/w )*( v(i)/w );
-			m_s = ( v(p) < 0 ? 1 : -1 ) * w * std::pow( m_s, 0.5 );
+			m_s = ( v(p) < 0 ? 1 : -1 ) * w * std::sqrt( m_s );
 		} else {
 			m_s = 0;
 		}
@@ -160,6 +161,19 @@
  *  @return \f$ h = v_p - s \f$
  */
 	inline const value_type h() const { return m_h; }
+/**
+ *  @return \f$ u \f$
+ */
+	inline vector_type u() const {
+		vector_type ret = zero_vector< value_type >( m_v.size() );
+		ret(m_p) = m_h;
+		project(ret,range(m_l,ret.size())) = project(m_v,range(m_l,m_v.size()));
+		return ret;
+	}
+/**
+ *  @return \f$ b \equiv \f$
+ */
+	inline const value_type b() const { return m_s * m_h; }
 };
 
 };

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/least_squares.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/least_squares.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/nnls.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/nnls.h Changed

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/qr_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/qr_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: qr_decomposition.h 66 2009-04-12 14:33:06Z matwey $
+    $Id: qr_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -19,31 +19,49 @@
 #ifndef _QR_DECOMPOSITION_H
 #define _QR_DECOMPOSITION_H
 
-#include <lsp/givens_rotation.h>
-
+#include <cmath>
 #include <limits>
 
 #include <boost/numeric/ublas/matrix.hpp>
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/storage.hpp>
 
-using namespace boost::numeric::ublas;
+#include <lsp/givens_rotation.h>
+#include <lsp/householder_transform.h>
 
 namespace lsp{
 
-namespace {
-	template<class T> static T shift( T q2, T q1, T e2, T e1 ) {
+using namespace boost::numeric::ublas;
+
+struct general_tag {};
+struct bidiag_tag {};
+struct tridiag_symmetric_tag {};
+
+template< class T, class Tag = general_tag > class qr_decomposition;
+
+namespace detail {
+	template<class T> T hypot( const T& x, const T& y ) {
 		typedef T value_type;
-		const value_type f = ( q2*q2 - q1*q1 +  e2*e2 - e1*e1 ) / ( 2*e2*q1 );
-		const value_type t = ( f < 0 ?
-              		- f + std::pow(( value_type(1)+f*f ),0.5) :
-              		- f - std::pow(( value_type(1)+f*f ),0.5) );
-			return ( q2*q2 + e2*(e2 - q1/t) );
-	};
+
+		if( std::abs(x) > std::abs(y) ) {
+			const value_type t = y/x;
+			return std::abs(x) * std::sqrt( value_type(1) + t*t );
+		}
+		if( y == value_type(0) ) {
+			return 0;
+		}
+		const value_type t = x/y;
+		return std::abs(y) * std::sqrt( value_type(1) + t*t );
+	}
+	template<class T> T wilkinson_shift( const T& a1, const T& b1, const T& a2 ) {
+		typedef T value_type;
+		const value_type d = (a1 - a2) / value_type(2);
+		const value_type q = hypot(d,b1);
+		return a2 - b1*b1 / ( d + ( d > 0 ? 1 : -1 ) * q );
+	}
 };
 
 /**
- *  @class qr_decomposition
  *  @brief A functor for the modified QR decomposition
  *
  *  Modified QR decomposition is iterative alghoritm that transformates
@@ -51,7 +69,7 @@
  *  Givens transformations. The decomposition is one of SVD parts.
  *
  */
-template<class T> class qr_decomposition {
+template<class T> class qr_decomposition<T, bidiag_tag> {
 public:
 	typedef T                                  matrix_type;   //!<
 	typedef typename matrix_type::value_type   value_type;    //!<
@@ -104,13 +122,18 @@
 	template<class M1, class M2> void apply( M1& left, M2& right, const range& cell, const regular_tag& ) const {
 		typedef givens_rotation< value_type > givens_rotation_type;
 
-		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_frobenius( m_matrix ) / m_matrix.size2();
+		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
 
 		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
-			while( std::abs( m_super( *it ) ) > lim ) {
-
-				value_type en1 = ( *it != cell(0) ? m_super(*it - 1) : value_type(0) );
-				value_type e0 = m_leading( cell(0) ) - shift( m_leading(*it + 1), m_leading(*it), m_super(*it), en1 ) / m_leading( cell(0) );
+			value_type last_super; /* use differential convergation control */
+			do {
+				last_super = m_super( *it );
+
+				const value_type en1 = ( *it != cell(0) ? m_super(*it - 1) : value_type(0) );
+				const value_type a1 = en1*en1 + m_leading(*it)*m_leading(*it);
+				const value_type a2 = m_super(*it)*m_super(*it) + m_leading(*it+1)*m_leading(*it+1);
+				const value_type b1 = m_super(*it)*m_leading(*it);
+				value_type e0 = m_leading( cell(0) ) - detail::wilkinson_shift( a1, b1, a2 ) / m_leading( cell(0) );
 				value_type z = m_super( cell(0) );
 
 				givens_rotation_type gr_left( e0, z );
@@ -132,7 +155,7 @@
 
 					gr_left = givens_rotation_type( m_super(*it2-1), z );
 				}
-			}
+			} while( std::abs( last_super - m_super( *it ) ) > lim );
 			m_super( *it ) = 0;
 		}
 	}
@@ -142,7 +165,7 @@
 		if( cell.size() == 1 ) /* Scalar is in diagonal form */
 			return ;
 
-		const value_type lim = 6 * std::numeric_limits< value_type >::epsilon() * norm_frobenius( m_matrix );
+		const value_type lim = 6 * std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
 
 		/* Looking for the zero diagonal element */
 		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend() ; ++it ) {
@@ -220,6 +243,240 @@
 
 };
 
+/**
+ *  @brief A functor for the modified QR decomposition
+ *
+ *  Modified QR decomposition is iterative alghoritm that transformates
+ *  tridiagonal banded symmetric matrix into diagonal form by the number of
+ *  Givens transformations. The decomposition is one of Schur decomposition parts.
+ *
+ */
+template<class T> class qr_decomposition<T, tridiag_symmetric_tag> {
+public:
+	typedef T                                  matrix_type;   //!<
+	typedef typename matrix_type::value_type   value_type;    //!<
+	typedef typename matrix_type::size_type    size_type;     //!<
+	typedef matrix_vector_slice< matrix_type > diagonal_type; //!<
+private:
+	struct regular_tag {};
+	struct left_tag {};
+	struct right_tag {};
+private:
+	matrix_type& m_matrix;
+	mutable diagonal_type m_super;
+	mutable diagonal_type m_leading;
+private:
+
+	template<class M1> void apply( M1& left, const range& cell, const regular_tag& ) const {
+		typedef givens_rotation< value_type > givens_rotation_type;
+
+		const value_type lim = std::numeric_limits< value_type >::epsilon() * norm_inf( m_matrix );
+
+		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
+			value_type last_super; /* use differential convergation control */
+			do {
+				last_super = m_super( *it );
+
+				value_type e0 = m_leading( cell(0) ) - detail::wilkinson_shift( m_leading(*it), m_super(*it), m_leading(*it+1) );
+				value_type z = m_super( cell(0) );
+
+				givens_rotation_type gr_left( e0, z );
+				for( range::const_iterator it2 = cell.begin() + 1; *it2 != *it + 2; ++it2 ) {
+
+					value_type t = m_super(*it2-1);
+					gr_left.apply( m_leading(*it2-1), t );
+					gr_left.apply( m_super(*it2-1), m_leading(*it2) );
+					gr_left.apply( m_leading(*it2-1), m_super(*it2-1) );
+					gr_left.apply( t, m_leading(*it2) );
+					gr_left.apply( row(left,*it2-1), row(left,*it2) );
+
+					if( *it2 == *it + 1 ) break;
+					z = gr_left.s() * m_super(*it2);
+					m_super(*it2) = gr_left.c() * m_super(*it2);
+
+					gr_left = givens_rotation_type( m_super(*it2-1), z );
+				}
+			} while( std::abs( last_super - m_super( *it ) ) > lim );
+			m_super( *it ) = 0;
+		}
+	}
+
+	template<class M1> void apply( M1& left, const range& cell ) const {
+
+		if( cell.size() == 1 ) /* Scalar is in diagonal form */
+			return ;
+
+		/* Looking for the zero upper diagonal element */
+		for( range::const_reverse_iterator it = cell.rbegin() + 1; it != cell.rend(); ++it ) {
+			if( m_super( *it ) == 0 ) {
+				apply( left, range( cell.start(), *it + 1 ) );
+				apply( left, range( *it + 1, cell.start() + cell.size() ) );
+				return ;
+			}
+		}
+
+		apply( left, cell, regular_tag() );
+
+	}
+

atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/schur_decomposition.h Added

@@ -0,0 +1,118 @@
+/*
+    $Id: schur_decomposition.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _SCHUR_DECOMPOSITION_H
+#define _SCHUR_DECOMPOSITION_H
+
+#include <boost/numeric/ublas/banded.hpp>
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+
+#include <lsp/tridiagonal_transform.h>
+#include <lsp/qr_decomposition.h>
+#include <lsp/utils.h>
+
+namespace lsp{
+
+using namespace boost::numeric::ublas;
+
+/**
+ *  @class schur_decomposition
+ *  @brief A functor for the Schur decomposition for symmetric matrices
+ *
+ *  Schur is a factorization of symmetric matrix that
+ *  \f
+ *  Q A Q^{T} = T \quad \mbox{where} \quad T = \left| \begin{array}{ccccc}
+ *     \lambda_1   &       &       &       &       \\
+ *           & \lambda_2   &       &       &       \\
+ *           &       &\ddots &       &       \\
+ *           &       &       &\lambda_{n-1}&       \\
+ *           &       &       &       & \lambda_n   \\
+ *  \end{array} \right|,\quad A \quad \mbox{is initial matrix,} \quad U \quad \mbox{is unitary matrix,} \quad \lambda_1 \ge \lambda_2 \ge \dots \ge \lambda_{n-1} \ge \lambda_{n}
+ *  \f
+ *
+ */
+template<class T> class schur_decomposition {
+public:
+	typedef T                                       matrix_type;               //!<
+	typedef typename matrix_type::value_type        value_type;                //!<
+	typedef typename matrix_type::size_type         size_type;                 //!<
+	typedef tridiagonal_transform< matrix_type >    tridiagonal_transform_type; //!<
+	typedef qr_decomposition< matrix_type, tridiag_symmetric_tag > qr_decomposition_type; //!<
+private:
+	matrix_type& m_matrix;
+	mutable tridiagonal_transform_type m_td_trans;
+	mutable qr_decomposition_type m_qr_decomp;
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix The reference to matrix object to be decomposited
+ *
+ *  Actual decomposition will be performed as soon as apply(M1& left) will be called.
+ *
+ */
+	schur_decomposition( matrix_type& matrix ):
+		m_matrix( matrix ),
+		m_td_trans( matrix ),
+		m_qr_decomp( matrix ) {
+	}
+
+/**
+ *  @brief Decomposition operaton
+ *  @paramout left  The left matrix
+ *
+ *  The routine decomposites the matrix.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+ *
+ *  \f$ M_{matrix} := T \f$
+ *
+ */
+	template<class M1> void apply( M1& left ) const {
+		typedef matrix_vector_slice< matrix_type > diagonal_type;
+		typedef vector< size_type > permutation_type;
+
+		m_td_trans.apply( left );
+
+		m_qr_decomp.apply( left );
+
+		diagonal_type lambda( m_matrix,
+			slice(0, 1, std::min( m_matrix.size1(), m_matrix.size2() )),
+			slice(0, 1, std::min( m_matrix.size1(), m_matrix.size2() )) );
+		permutation_type pm( lambda.size() );
+		for( typename permutation_type::iterator it = pm.begin(); it != pm.end(); ++it ){
+			*it = it.index();
+		}
+
+		std::sort( pm.begin(), pm.end(), vector_less< diagonal_type, std::greater< typename diagonal_type::value_type > >( lambda ) );
+
+		for( typename permutation_type::size_type it = 0; it != pm.size(); ++it ){
+			if( it < pm(it) ) {
+				row(m_matrix, pm(it)).swap( row(m_matrix, it) );
+				row(left, pm(it)).swap( row(left, it) );
+			}
+		}
+	}
+};
+
+
+};
+
+#endif // _SCHUR_DECOMPOSITION_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/singular_decomposition.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/singular_decomposition.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: singular_decomposition.h 55 2009-02-05 11:01:48Z matwey $
+    $Id: singular_decomposition.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -28,10 +28,10 @@
 #include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/numeric/ublas/banded.hpp>
 
-using namespace boost::numeric::ublas;
-
 namespace lsp{
 
+using namespace boost::numeric::ublas;
+
 /**
  *  @class singular_decomposition
  *  @brief A functor for the singular value decomposition (SVD)
@@ -57,7 +57,7 @@
 	typedef typename matrix_type::size_type         size_type;                 //!< The type for seeking in the matrix object
 	typedef banded_adaptor< matrix_type >           banded_adaptor_type;       //!< The type of banded adaptor
 	typedef bidiagonal_transform< matrix_type >     bidiagonal_transform_type; //!< The type of a functor for bidiagonal transformation
-	typedef qr_decomposition< banded_adaptor_type > qr_decomposition_type;     //!< The type of a functor for QR decomposition
+	typedef qr_decomposition< banded_adaptor_type, bidiag_tag > qr_decomposition_type;     //!< The type of a functor for QR decomposition
 private:
 	matrix_type& m_matrix;
 	mutable bidiagonal_transform_type m_bd_trans;

atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/tridiagonal_transform.h Added

@@ -0,0 +1,115 @@
+/*
+    $Id: tridiagonal_transform.h 80 2010-12-10 17:11:54Z matwey $
+    Copyright (C) 2010  Matwey V. Kornilov <matwey.kornilov@gmail.com>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef _TRIDIAGONAL_TRANSFORM_H
+#define _TRIDIAGONAL_TRANSFORM_H
+
+#include <limits>
+
+#include <boost/numeric/ublas/matrix.hpp>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
+#include <boost/numeric/ublas/vector.hpp>
+#include <boost/numeric/ublas/vector_proxy.hpp>
+
+#include <lsp/householder_transform.h>
+
+namespace lsp{
+
+using namespace boost::numeric::ublas;
+
+namespace detail {
+	template<class M> symmetric_adaptor<M> make_symmetric_adaptor( M& e) {
+		return symmetric_adaptor<M>(e);
+	}
+	template<class M> symmetric_adaptor<const M> make_symmetric_adaptor( const M& e) {
+		return symmetric_adaptor<const M>(e);
+	}
+}
+	
+/**
+ *  @class trdiagonal_transform
+ *  @brief A functor for the transformation symmetric matrix into the tridiagonal form using Householder transformations.
+ *
+ */
+template<class T> class tridiagonal_transform {
+public:
+	typedef T                                matrix_type; //!< The type of the matrix object to be trasformed
+	typedef typename matrix_type::value_type value_type;  //!< The type of the elements stored in the matrix_type
+	typedef typename matrix_type::size_type  size_type;   //!< The type for seeking in the matrix object
+
+private:
+	matrix_type& m_matrix;
+	size_type    m_size;
+
+public:
+/**
+ *  @brief An object constructor
+ *  @paramin,out matrix The reference to matrix object to be transformed
+ *
+ *  Actual transformation will be performed as soon as apply(M1& left, M2& right) will be called.
+ *
+ */
+	tridiagonal_transform( matrix_type& matrix ):
+		m_matrix( matrix ),
+		m_size( matrix.size1() ) {
+
+		assert( matrix.size1() == matrix.size2() );
+	}
+
+/**
+ *  @brief Transformation operaton
+ *  @paramout left  The left matrix
+ *
+ *  The routine calculates and makes transformation.
+ *  Intrinsic assumption is that the all matrix are size-suitable.
+ *
+ *  \f$ M_{left} := Q M_{left} \f$
+  *
+ *  \f$ M_{matrix} := Q M_{matrix} Q^{T} \equiv B \f$
+ *
+ */
+	template<class M1> void apply( M1& left ) const {
+		typedef vector< value_type >                 vector_type;
+		typedef householder_transform< vector_type > householder_transform_type;
+
+		assert( left.size1()  == m_matrix.size1() );
+
+		if( m_size < 3 )
+			return ;
+
+		for( size_type i = 0; i < m_size - 1; ++i ) {
+
+			householder_transform_type hleft( i+2, i+1, column(m_matrix,i) );
+			hleft.apply( column(m_matrix,i) );
+			hleft.apply( left, row_major_tag() );
+			
+			matrix_range< matrix_type > sub_submatrix( m_matrix, range(i+1, m_size), range(i+1, m_size) );
+			value_type  c = value_type(1) / hleft.b();
+			vector_type u = project( hleft.u(), range(i+1,m_size) );
+			vector_type p = prod( sub_submatrix, u );
+			vector_type w = p + c / value_type(2) * inner_prod(p,u) * u;
+
+			sub_submatrix += detail::make_symmetric_adaptor( c * (outer_prod(u,w) + outer_prod(w,u)) );
+		}
+	}
+
+};
+
+};
+
+#endif // _TRIDIAGONAL_TRANSFORM_H

atmos-2.96.9.tar.gz/third-party/lsp/include/lsp/utils.h -> atmos-2.97.1.tar.gz/third-party/lsp/include/lsp/utils.h Changed

@@ -1,5 +1,5 @@
 /*
-    $Id: utils.h 59 2009-03-01 18:25:07Z matwey $
+    $Id: utils.h 80 2010-12-10 17:11:54Z matwey $
     Copyright (C) 2008  Matwey V. Kornilov <matwey.kornilov@gmail.com>
 
     This program is free software: you can redistribute it and/or modify
@@ -24,10 +24,17 @@
 #include <cmath>
 #include <functional>
 
+#include <boost/numeric/ublas/expression_types.hpp>
+#include <boost/type_traits.hpp>
+#include <boost/utility/enable_if.hpp>
+
 /**
  *  @brief Basic namespace
  */
 namespace lsp {
+
+using namespace boost::numeric::ublas;
+
 /**
  *  @class less_abs
  *  @brief Comparsion of absoulte values
@@ -170,6 +177,19 @@
 
 };
 
+template<class T, class Enable = void> struct temporary_type_traits {
+	typedef T type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, scalar_tag > >::type > {
+	typedef typename T::expression_type::value_type type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, vector_tag > >::type > {
+	typedef typename T::vector_temporary_type type;
+};
+template<class T> struct temporary_type_traits<T, typename boost::enable_if< boost::is_same< typename T::type_category, matrix_tag > >::type > {
+	typedef typename T::matrix_temporary_type type;
+};
+
 };
 
 #endif // _UTILS_H