Changeset 1715 for trunk/src/STAtmosphere.cpp
- Timestamp:
- 03/17/10 22:03:25 (14 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/STAtmosphere.cpp
r1713 r1715 65 65 **/ 66 66 STAtmosphere::STAtmosphere(double wvScale, double maxAlt, size_t nLayers) : 67 itsGndTemperature(288.), itsGndPressure(101325.), itsGndHumidity(0.5),68 itsLapseRate(0.0065), itsWVScale(wvScale), itsMaxAlt(maxAlt),69 67 itsHeights(nLayers), itsTemperatures(nLayers), 70 itsDryPressures(nLayers), itsVapourPressures(nLayers) 68 itsDryPressures(nLayers), itsVapourPressures(nLayers), 69 itsGndTemperature(288.), itsPressure(101325.), itsGndHumidity(0.5), 70 itsLapseRate(0.0065), itsWVScale(wvScale), itsMaxAlt(maxAlt), itsObsHeight(200.) 71 71 { 72 72 recomputeAtmosphereModel(); … … 76 76 * Constructor with explicitly given parameters of the atmosphere 77 77 * @param[in] temperature air temperature at the observatory (K) 78 * @param[in] pressure air pressure at the observatory(Pascals)78 * @param[in] pressure air pressure at the sea level (Pascals) 79 79 * @param[in] humidity air humidity at the observatory (fraction) 80 80 * @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, 88 88 double wvScale, double maxAlt, size_t nLayers) : 89 itsGndTemperature(temperature), itsGndPressure(pressure), itsGndHumidity(humidity),90 itsLapseRate(lapseRate), itsWVScale(wvScale), itsMaxAlt(maxAlt),91 89 itsHeights(nLayers), itsTemperatures(nLayers), 92 itsDryPressures(nLayers), itsVapourPressures(nLayers) 90 itsDryPressures(nLayers), itsVapourPressures(nLayers), 91 itsGndTemperature(temperature), itsPressure(pressure), itsGndHumidity(humidity), 92 itsLapseRate(lapseRate), itsWVScale(wvScale), itsMaxAlt(maxAlt), itsObsHeight(200.) 93 93 { 94 94 recomputeAtmosphereModel(); … … 98 98 * Set the new weather station data, recompute the model 99 99 * @param[in] temperature air temperature at the observatory (K) 100 * @param[in] pressure air pressure at the observatory(Pascals)100 * @param[in] pressure air pressure at the sea level (Pascals) 101 101 * @param[in] humidity air humidity at the observatory (fraction) 102 102 **/ … … 104 104 { 105 105 itsGndTemperature = temperature; 106 its GndPressure = pressure;106 itsPressure = pressure; 107 107 itsGndHumidity = humidity; 108 108 recomputeAtmosphereModel(); … … 110 110 111 111 /** 112 * Set the elevation of the observatory (height above mean sea level) 113 * By default, 200m is assumed. 114 * @param[in] elev elevation in metres 115 **/ 116 void STAtmosphere::setObservatoryElevation(double elev) 117 { 118 itsObsHeight = elev; 119 recomputeAtmosphereModel(); 120 } 121 122 123 /** 112 124 * Build the atmosphere model based on exponential fall-off, ideal gas and hydrostatic 113 125 * equilibrium. The model parameters are taken from the data members of this class. … … 116 128 { 117 129 AlwaysAssert(itsGndTemperature > 0, AipsError); 118 AlwaysAssert(its GndPressure > 0., AipsError);130 AlwaysAssert(itsPressure > 0., AipsError); 119 131 AlwaysAssert((itsGndHumidity >= 0.) && (itsGndHumidity<=1.), AipsError); 120 132 AlwaysAssert(itsMaxAlt > 0., AipsError); … … 126 138 // free-fall acceleration 127 139 const double g = 9.81; 140 128 141 const double wvGndSaturationPressure = wvSaturationPressure(itsGndTemperature); 142 const double gndPressure = itsPressure*exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)* 143 (itsObsHeight+0.5*itsLapseRate*itsObsHeight*itsObsHeight/itsGndTemperature)); 129 144 for (size_t layer = 0; layer < nLayers(); ++layer) { 130 145 const double height = double(layer)*heightStep; 131 146 itsHeights[layer] = height; 132 147 itsTemperatures[layer] = itsGndTemperature/(1.+itsLapseRate*height/itsGndTemperature); 133 const double pressure = itsGndPressure * exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)*148 const double pressure = gndPressure * exp(-M*g/(QC::R.get().getValue()*itsGndTemperature)* 134 149 (height+0.5*itsLapseRate*height*height/itsGndTemperature)); 135 150 itsVapourPressures[layer] = casa::min(itsGndHumidity*exp(-height/itsWVScale)*wvGndSaturationPressure, 136 151 wvSaturationPressure(itsTemperatures[layer])); 137 152 itsDryPressures[layer] = pressure - itsVapourPressures[layer]; 153 std::cout<<"layer="<<layer<<": H="<<itsHeights[layer]<<" T="<<itsTemperatures[layer]<< 154 " Pvap="<<itsVapourPressures[layer]<<" Pdry="<<itsDryPressures[layer]<<endl; 138 155 } 139 156 } … … 166 183 double STAtmosphere::wvSaturationPressure(double temperature) 167 184 { 168 if (temperature >215.) {185 if (temperature <= 215.) { 169 186 return 0.; 170 187 }
Note: See TracChangeset
for help on using the changeset viewer.