Line data    Source code

       1             :  /*******************************************************************************
       2             :  * ALMA - Atacama Large Millimeter Array
       3             :  * (c) Instituto de Estructura de la Materia, 2011
       4             :  * (in the framework of the ALMA collaboration).
       5             :  * All rights reserved.
       6             :  * 
       7             :  * This library is free software; you can redistribute it and/or
       8             :  * modify it under the terms of the GNU Lesser General Public
       9             :  * License as published by the Free Software Foundation; either
      10             :  * version 2.1 of the License, or (at your option) any later version.
      11             :  * 
      12             :  * This library is distributed in the hope that it will be useful, 
      13             :  * but WITHOUT ANY WARRANTY; without even the implied warranty of
      14             :  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
      15             :  * Lesser General Public License for more details.
      16             :  * 
      17             :  * You should have received a copy of the GNU Lesser General Public
      18             :  * License along with this library; if not, write to the Free Software
      19             :  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA
      20             :  *******************************************************************************/
      21             : 
      22             : #include <string>
      23             : #include <vector>
      24             : #include <iostream>
      25             : #include <fstream>
      26             : #include <math.h> 
      27             : #include <string.h>
      28             : #include <stdlib.h>
      29             : #include <iomanip>
      30             : 
      31             : using namespace std;
      32             : 
      33             : #include <atmosphere/ATM/ATMPercent.h>
      34             : #include <atmosphere/ATM/ATMPressure.h>
      35             : #include <atmosphere/ATM/ATMNumberDensity.h>
      36             : #include <atmosphere/ATM/ATMMassDensity.h>
      37             : #include <atmosphere/ATM/ATMTemperature.h>
      38             : #include <atmosphere/ATM/ATMLength.h>
      39             : #include <atmosphere/ATM/ATMInverseLength.h>
      40             : #include <atmosphere/ATM/ATMOpacity.h>
      41             : #include <atmosphere/ATM/ATMAngle.h>
      42             : #include <atmosphere/ATM/ATMHumidity.h>
      43             : #include <atmosphere/ATM/ATMFrequency.h>
      44             : #include <atmosphere/ATM/ATMWaterVaporRadiometer.h>
      45             : #include <atmosphere/ATM/ATMWVRMeasurement.h>
      46             : #include <atmosphere/ATM/ATMProfile.h>
      47             : #include <atmosphere/ATM/ATMSpectralGrid.h>
      48             : #include <atmosphere/ATM/ATMRefractiveIndex.h>
      49             : #include <atmosphere/ATM/ATMRefractiveIndexProfile.h>
      50             : #include <atmosphere/ATM/ATMSkyStatus.h>
      51             : using namespace atm;
      52             : 
      53           1 : int main()
      54             : {  
      55           1 :   unsigned int atmType = 1; // TROPICAL  // Atmospheric type (to reproduce behavior above the tropopause)
      56           2 :   Temperature      T( 270.0,"K" );       // Ground temperature
      57           2 :   Pressure         P( 560.0,"mb");       // Ground Pressure
      58           2 :   Humidity         H(  20.0,"%" );       // Ground Relative Humidity (indication)
      59           2 :   Length         Alt(  5000,"m" );       // Altitude of the site 
      60           2 :   Length         WVL(   2.0,"km");       // Water vapor scale height
      61           1 :   double         TLR=  -5.6      ;       // Tropospheric lapse rate (must be in K/km)
      62           2 :   Length      topAtm(  95.0,"km");       // Upper atm. boundary for calculations
      63           2 :   Pressure     Pstep(   5.0,"mb");       // Primary pressure step (10.0 mb)
      64             : 
      65           1 :   string namefile;
      66           1 :   namefile="results.txt" ;
      67           1 :   ofstream myfile(namefile.c_str());
      68             : 
      69             : 
      70           1 :   myfile<<"! SkyStatusTest: BASIC ATMOSPHERIC PARAMETERS TO GENERATE REFERENCE ATMOSPHERIC PROFILE"<<endl; 
      71           1 :   myfile<<"! SkyStatusTest:   "<<endl;
      72           1 :   myfile<<"! SkyStatusTest: Ground temperature T:         " << T.get()         << " K"    <<endl; 
      73           1 :   myfile<<"! SkyStatusTest: Ground pressure P:            " << P.get("mb")     << " mb"   <<endl; 
      74           1 :   myfile<<"! SkyStatusTest: Relative humidity rh:         " << H.get("%")      << " %"    <<endl; 
      75           1 :   myfile<<"! SkyStatusTest: Scale height h0:              " << WVL.get("km")   << " km"   <<endl; 
      76           1 :   myfile<<"! SkyStatusTest: Pressure step dp:             " << Pstep.get("mb") << " mb"   <<endl; 
      77           1 :   myfile<<"! SkyStatusTest: Altitude alti:                " << Alt.get()       << " m"    <<endl; 
      78           1 :   myfile<<"! SkyStatusTest: Attitude top atm profile atmh:" << topAtm.get("km")<< " km"   <<endl; 
      79           1 :   myfile<<"! SkyStatusTest: Tropospherique lapse rate:    " << TLR                << " K/km" <<endl;
      80             : 
      81           1 :   double pressureStepFactor = 1.075;
      82           1 :   AtmProfile myProfile( Alt, P, T, TLR, H, WVL, Pstep, pressureStepFactor,  topAtm, atmType );
      83             : 
      84           1 :   myfile<<"! SkyStatusTest: Object myProfile built with the AtmProfile CONSTRUCTOR and the above entries"<< endl; 
      85           1 :   myfile<<"! SkyStatusTest: Atmospheric type:             " << myProfile.getAtmosphereType() <<endl;
      86           1 :   myfile<<"! SkyStatusTest: Number of layers returned:  " << myProfile.getNumLayer() <<endl;
      87           1 :   myfile<<"! SkyStatusTest: Layer parameters:  " <<endl;
      88             :    
      89             :     
      90           1 :     myfile << "! SkyStatusTest: size of Temperature Profile: " << myProfile.getTemperatureProfile().size() << endl;
      91           1 :     myfile << "!   H0_layer  H1_layer  Layer Thickness   Pressure    Temperature             WaterVapor             O3           CO          N2O          HCl          HCN  "  << endl;
      92           1 :     myfile << "!      (m)       (m)          (m)           (mb)          (K)            kg m^-3        m^-3        m^-3         m^-3         m^-3         m^-3         m^-3 "  << endl;
      93           1 :     double alt0=Alt.get();
      94           1 :     double alt1=alt0;
      95          43 :    for(unsigned int i=0; i<myProfile.getNumLayer(); i++){
      96          42 :      alt1=alt0+myProfile.getLayerThickness(i).get("m");
      97          42 :      myfile << "!"
      98          42 :             << fixed << showpoint << setprecision(3) << setw(11) << alt0 
      99          42 :             << fixed << showpoint << setprecision(3) << setw(11) << alt1
     100          84 :             << fixed << showpoint << setprecision(3) << setw(11) << myProfile.getLayerThickness(i).get("m")  << "      "       
     101          84 :             << fixed << showpoint << setprecision(3) << setw(9) << myProfile.getLayerPressure(i).get("mb")  << "    "
     102          84 :             << fixed << showpoint << setprecision(3) << setw(9) << myProfile.getLayerTemperature(i).get("K") << "      "  
     103          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerWaterVaporMassDensity(i).get("kgm**-3")  
     104          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerWaterVaporNumberDensity(i).get("m**-3")  
     105          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerO3(i).get("m**-3")    
     106          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerCO(i).get("m**-3")        
     107          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerN2O(i).get("m**-3")
     108          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerHCl(i).get("m**-3")
     109          84 :             << scientific << showpoint << setprecision(5) << setw(13) << myProfile.getLayerHCN(i).get("m**-3")
     110          42 :             << endl;
     111          42 :      alt0=alt1;
     112             :    }
     113             : 
     114           1 :   myfile   << "! SkyStatusTest: First guess precipitable water vapor content: " << myProfile.getGroundWH2O().get("mm") << "mm" << endl;
     115           1 :   myfile   << "! SkyStatusTest: (This value is estimated from the relative humidity at ground level and the water vapor scale height)" << endl;
     116           1 :   myfile   << "! SkyStatusTest:   " << endl;
     117             :     
     118           1 :   myfile   << "! SkyStatusTest:    TEST FOR ONE SPECTRAL WINDOW WITH SEVERAL CHANNELS     " << endl;
     119           1 :   myfile   << "! SkyStatusTest:    =====================================================  " << endl;
     120           1 :   myfile   << "! SkyStatusTest:  " << endl;
     121             : 
     122           1 :   unsigned int numchan=100;
     123           1 :   unsigned int refchan=51;      
     124           2 :   Frequency reffreq(500.10,"GHz");  
     125             :   /* unsigned int numchan=50;
     126             :   unsigned int refchan=26;      
     127             :   Frequency reffreq(372.00,"GHz"); */
     128           2 :   Frequency chansep(0.1,"GHz");
     129             :   
     130           1 :   SpectralGrid spectral_band(numchan, refchan, reffreq, chansep);
     131           1 :   RefractiveIndexProfile refractive_index_spectrum(spectral_band, myProfile); 
     132           1 :   SkyStatus sky_spectrum(refractive_index_spectrum);
     133             : 
     134           1 :   sky_spectrum.setAirMass(1.0);
     135           1 :   sky_spectrum.setUserWH2O(0.33,"mm");
     136             : 
     137           1 :   myfile << "! SkyStatusTest: water vapor column: " << sky_spectrum.getUserWH2O().get("mm") << " mm" << endl;
     138           1 :   myfile << "! SkyStatusTest: Air mass          : " << sky_spectrum.getAirMass() << endl;
     139           1 :   myfile << "! " << endl;
     140             : 
     141             : 
     142             : 
     143             : 
     144             :   /*  myfile << "! Frequency         Sky        D_Sky_Trm     T_EBB     < --------------------------------------------------------------------------------------" <<
     145             :     " Total opacities along the atmospheric path ------------------------------------------------------------------------------" <<
     146             :     "----------------------------------------------------------------------->" << endl;
     147             :   myfile << "!   (GHz)       Transmission  (D_nu=100Hz)    (K)       Wet_CIA      Dry_CIA      H2O Line      HH16O_Lin      HH16OV2Lin        HH17O_Lin         HH18O_Lin  " <<
     148             :     "        HDO_Lin           O2_Lines           16O16O_Lin           16O2Vib_L              16O17O_Lin             16O18O_Lin                 O3_Lines   16O3_Lin   16O3_V1_L " <<
     149             :     " 16O3_V2_L  16O3_V3_L  16O16O17O  16O16O18O  16O17O16O  16O18O16O     N2O        CO        HCl        HCN" << endl;
     150             :   */
     151             :   
     152           1 :   myfile << "! Frequency         Sky        D_Sky_Trm     T_EBB     < -------------------------------------------------------------------------------------- Total opacities along the atmospheric path ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------->" << endl;
     153           1 :   myfile << "!   (GHz)       Transmission  (D_nu=100Hz)    (K)       Wet_CIA    Dry_CIA     H2O Line    HH16O_Lin  HH16OV2Lin  HH17O_Lin    HH18O_Lin    HDO_Lin    O2_Lines   16O16O_Lin   16O2Vib_L  16O17O_Lin  16O18O_Lin   O3_Lines    16O3_Lin   16O3_V1_L   16O3_V2_L   16O3_V3_L   16O16O17O   16O16O18O   16O17O16O   16O18O16O       N2O         CO          HCl         HCN  " << endl;
     154             : 
     155           1 :   unsigned int numchanDF=3;
     156           1 :   unsigned int refchanDF=2;        
     157           2 :   Frequency chansepDF(0.0000005,"GHz");
     158           1 :   Frequency reffreqDF;
     159             :   double TransmissionVariation;
     160             :   double delta_right;
     161             :   double delta_left;
     162             : 
     163             : 
     164          99 :   for(unsigned int i=1; i<sky_spectrum.getNumChan(0)-1; i++){
     165             : 
     166          98 :     reffreqDF=sky_spectrum.getChanFreq(i);
     167          98 :     SpectralGrid spectral_bandDF(numchanDF, refchanDF, reffreqDF, chansepDF);
     168          98 :     RefractiveIndexProfile refractive_index_spectrumDF(spectral_bandDF, myProfile); 
     169          98 :     SkyStatus sky_spectrumDF(refractive_index_spectrumDF);
     170          98 :     delta_right=(sky_spectrumDF.getSkyTransmission(2)-sky_spectrumDF.getSkyTransmission(1)); // /sky_spectrumDF.getSkyTransmission(1);
     171          98 :     delta_left= (sky_spectrumDF.getSkyTransmission(0)-sky_spectrumDF.getSkyTransmission(1)); // /sky_spectrumDF.getSkyTransmission(1);;
     172          98 :     if(abs(delta_right)>abs(delta_left)){TransmissionVariation=delta_right;}else{TransmissionVariation=delta_left;}   
     173             :     
     174             :     //    if(int(10*sky_spectrum.getNumChan(0)/i)==10*int(sky_spectrum.getNumChan(0)/i)){cout << 10*int(sky_spectrum.getNumChan(0)/i) << "% completed" << endl;}
     175             :     
     176         196 :     myfile << fixed << showpoint << setprecision(7) << setw(13) << sky_spectrum.getChanFreq(i).get("GHz") 
     177          98 :            << fixed << showpoint << setprecision(8) << setw(14) << sky_spectrum.getSkyTransmission(i)
     178          98 :            << fixed << showpoint << setprecision(8) << setw(14) << TransmissionVariation
     179         196 :            << fixed << showpoint << setprecision(4) << setw(10) << sky_spectrum.getTebbSky(i).get("K") 
     180         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getH2OContOpacity(i).get()  
     181         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getDryContOpacity(i).get() 
     182         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getH2OLinesOpacity(i).get() 
     183         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHH16OLinesOpacity(i).get() 
     184         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHH16OV2LinesOpacity(i).get() 
     185         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHH17OLinesOpacity(i).get() 
     186         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHH18OLinesOpacity(i).get() 
     187         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHDOLinesOpacity(i).get() 
     188         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getO2LinesOpacity(i).get() 
     189         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16OLinesOpacity(i).get() 
     190         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16OVIBLinesOpacity(i).get() 
     191         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O17OLinesOpacity(i).get() 
     192         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O18OLinesOpacity(i).get() 
     193         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getO3LinesOpacity(i).get() 
     194         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O16OLinesOpacity(i).get() 
     195         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O16OV1LinesOpacity(i).get() 
     196         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O16OV2LinesOpacity(i).get() 
     197         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O16OV3LinesOpacity(i).get() 
     198         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O17OLinesOpacity(i).get() 
     199         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O16O18OLinesOpacity(i).get() 
     200         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O17O16OLinesOpacity(i).get() 
     201         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.get16O18O16OLinesOpacity(i).get() 
     202         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getN2OLinesOpacity(i).get() 
     203         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getCOLinesOpacity(i).get() 
     204         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHClLinesOpacity(i).get() 
     205         196 :            << fixed << showpoint << setprecision(6) << setw(12) << sky_spectrum.getHCNLinesOpacity(i).get() 
     206          98 :            << endl;    
     207          98 :   }
     208             : 
     209           1 :   return 0;
     210             : 
     211           1 : }
     212             : 
     213             : 
     214             : 
     215             : 
     216             : 
     217             : 
     218             :