Changeset 1715 for trunk

Timestamp:

03/17/10 22:03:25 (15 years ago)

Author:

Max Voronkov

Message:

changed pressure to the mean sea level to match miriad, bugfix, C++ code is now exposed to python

Location:

trunk/src

Files:

• STAtmosphere.cpp (modified) (9 diffs)
• STAtmosphere.h (modified) (4 diffs)
• python_STAtmosphere.cpp (added)
• python_asap.cpp (modified) (1 diff)
• python_asap.h (modified) (1 diff)

trunk/src/STAtmosphere.cpp

@@ -65 +65 @@
  **/
 STAtmosphere::STAtmosphere(double wvScale, double maxAlt, size_t nLayers) :
-   itsGndTemperature(288.), itsGndPressure(101325.), itsGndHumidity(0.5), 
-   itsLapseRate(0.0065), itsWVScale(wvScale), itsMaxAlt(maxAlt),
    itsHeights(nLayers), itsTemperatures(nLayers), 
-   itsDryPressures(nLayers), itsVapourPressures(nLayers) 
+   itsDryPressures(nLayers), itsVapourPressures(nLayers),
+   itsGndTemperature(288.), itsPressure(101325.), itsGndHumidity(0.5), 
+   itsLapseRate(0.0065), itsWVScale(wvScale), itsMaxAlt(maxAlt), itsObsHeight(200.)
 {
   recomputeAtmosphereModel();
@@ -76 +76 @@
  * Constructor with explicitly given parameters of the atmosphere 
  * @param[in] temperature air temperature at the observatory (K)
- * @param[in] pressure air pressure at the observatory (Pascals)
+ * @param[in] pressure air pressure at the sea level (Pascals)
  * @param[in] humidity air humidity at the observatory (fraction)
  * @param[in] lapseRate temperature lapse rate (K/m), default is 0.0065 K/m to match MIRIAD and ISA 
@@ -87 +87 @@
 STAtmosphere::STAtmosphere(double temperature, double pressure, double humidity, double lapseRate, 
                double wvScale, double maxAlt, size_t nLayers) :
-   itsGndTemperature(temperature), itsGndPressure(pressure), itsGndHumidity(humidity), 
-   itsLapseRate(lapseRate), itsWVScale(wvScale), itsMaxAlt(maxAlt),
    itsHeights(nLayers), itsTemperatures(nLayers), 
-   itsDryPressures(nLayers), itsVapourPressures(nLayers) 
+   itsDryPressures(nLayers), itsVapourPressures(nLayers),
+   itsGndTemperature(temperature), itsPressure(pressure), itsGndHumidity(humidity), 
+   itsLapseRate(lapseRate), itsWVScale(wvScale), itsMaxAlt(maxAlt), itsObsHeight(200.)
 {
   recomputeAtmosphereModel();
@@ -98 +98 @@
  * Set the new weather station data, recompute the model 
  * @param[in] temperature air temperature at the observatory (K)
- * @param[in] pressure air pressure at the observatory (Pascals)
+ * @param[in] pressure air pressure at the sea level (Pascals)
  * @param[in] humidity air humidity at the observatory (fraction)
  **/
@@ -104 +104 @@
 {
   itsGndTemperature = temperature;
-  itsGndPressure = pressure;
+  itsPressure = pressure;
   itsGndHumidity = humidity;
   recomputeAtmosphereModel();
@@ -110 +110 @@
 
 /**
+ * Set the elevation of the observatory (height above mean sea level)
+ * By default, 200m is assumed.
+ * @param[in] elev elevation in metres
+ **/
+void STAtmosphere::setObservatoryElevation(double elev)
+{
+  itsObsHeight = elev;
+  recomputeAtmosphereModel();  
+}
+
+
+/**
  * Build the atmosphere model based on exponential fall-off, ideal gas and hydrostatic
  * equilibrium. The model parameters are taken from the data members of this class.
@@ -116 +128 @@
 {
   AlwaysAssert(itsGndTemperature > 0, AipsError);
-  AlwaysAssert(itsGndPressure > 0., AipsError);
+  AlwaysAssert(itsPressure > 0., AipsError);
   AlwaysAssert((itsGndHumidity >= 0.) && (itsGndHumidity<=1.), AipsError);
   AlwaysAssert(itsMaxAlt > 0., AipsError);
@@ -126 +138 @@
   // free-fall acceleration
   const double g = 9.81;
+  
   const double wvGndSaturationPressure = wvSaturationPressure(itsGndTemperature);
+  const double gndPressure = itsPressure*exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)*
+                   (itsObsHeight+0.5*itsLapseRate*itsObsHeight*itsObsHeight/itsGndTemperature));
   for (size_t layer = 0; layer < nLayers(); ++layer) {
        const double height = double(layer)*heightStep;
        itsHeights[layer] = height;
        itsTemperatures[layer] = itsGndTemperature/(1.+itsLapseRate*height/itsGndTemperature);
-       const double pressure = itsGndPressure * exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)*
+       const double pressure = gndPressure * exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)*
                    (height+0.5*itsLapseRate*height*height/itsGndTemperature));
        itsVapourPressures[layer] = casa::min(itsGndHumidity*exp(-height/itsWVScale)*wvGndSaturationPressure,
                                              wvSaturationPressure(itsTemperatures[layer]));
        itsDryPressures[layer] = pressure - itsVapourPressures[layer];                                      
+       std::cout<<"layer="<<layer<<": H="<<itsHeights[layer]<<" T="<<itsTemperatures[layer]<<
+           " Pvap="<<itsVapourPressures[layer]<<" Pdry="<<itsDryPressures[layer]<<endl;
   }
 }
@@ -166 +183 @@
 double STAtmosphere::wvSaturationPressure(double temperature)
 {
-  if (temperature > 215.) {
+  if (temperature <= 215.) {
       return 0.;
   }

trunk/src/STAtmosphere.h

@@ -110 +110 @@
 
   /**
+   * Set the elevation of the observatory (height above mean sea level)
+   * By default, 200m is assumed.
+   * @param[in] elev elevation in metres
+   **/
+  void setObservatoryElevation(double elev);
+
+  /**
    * Calculate zenith opacity at the given frequency. This is a simplified version
    * of the routine implemented in MIRIAD, which calculates just zenith opacity and
@@ -136 +143 @@
    * @return zenith opacity (nepers, i.e. dimensionless)
    **/ 
-   double opacity(double freq, double el) const;
+  double opacity(double freq, double el) const;
 
   /**
@@ -226 +233 @@
   double itsGndTemperature;
   
-  // ground level pressure (Pascals)
-  double itsGndPressure;
+  // sea level pressure (Pascals)
+  double itsPressure;
   
   // ground level humidity (fraction)
@@ -240 +247 @@
   // altitude of the highest layer of the model (m)
   double itsMaxAlt;
+  
+  // observatory elevation (m)
+  double itsObsHeight;
 };
 

trunk/src/python_asap.cpp

@@ -76 +76 @@
   asap::python::python_Logger();
   asap::python::python_STCoordinate();
+  asap::python::python_STAtmosphere();
 
 #ifndef HAVE_LIBPYRAP

trunk/src/python_asap.h

@@ -48 +48 @@
     void python_Logger();
     void python_STCoordinate();
+    void python_STAtmosphere();
 
   } // python