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 <sstream>
      27             : #include <math.h>
      28             : 
      29             : using namespace std;
      30             : 
      31             : #include <atmosphere/ATM/ATMPercent.h>
      32             : #include <atmosphere/ATM/ATMPressure.h>
      33             : #include <atmosphere/ATM/ATMNumberDensity.h>
      34             : #include <atmosphere/ATM/ATMMassDensity.h>
      35             : #include <atmosphere/ATM/ATMTemperature.h>
      36             : #include <atmosphere/ATM/ATMLength.h>
      37             : #include <atmosphere/ATM/ATMInverseLength.h>
      38             : #include <atmosphere/ATM/ATMOpacity.h>
      39             : #include <atmosphere/ATM/ATMAngle.h>
      40             : #include <atmosphere/ATM/ATMHumidity.h>
      41             : #include <atmosphere/ATM/ATMFrequency.h>
      42             : #include <atmosphere/ATM/ATMWaterVaporRadiometer.h>
      43             : #include <atmosphere/ATM/ATMWVRMeasurement.h>
      44             : #include <atmosphere/ATM/ATMProfile.h>
      45             : #include <atmosphere/ATM/ATMSpectralGrid.h>
      46             : #include <atmosphere/ATM/ATMRefractiveIndex.h>
      47             : #include <atmosphere/ATM/ATMRefractiveIndexProfile.h>
      48             : #include <atmosphere/ATM/ATMSkyStatus.h>
      49             : 
      50           1 : int main()
      51             : {  
      52             :   // Initialize the Atmospheric model
      53             : 
      54           1 :   unsigned int      atmType = 1;       // 1=tropical (ALMA site), 2=midlatSummer, 3=midlatWinter 
      55           3 :   atm::Temperature  T(273. ,"K" );     // Ground temperature
      56           3 :   atm::Pressure     P(55000. ,"Pa");   // Ground Pressure
      57           3 :   atm::Humidity     H(8.,"%" );       // Ground Relative Humidity (indication)
      58           3 :   atm::Length       Alt(5000,"m" );    // Altitude of the site
      59           3 :   atm::Length       WVL(   3.0,"km");  // Water vapor scale height
      60           1 :   double            TLR=  -6.5      ;  // Tropospheric lapse rate (must be in K/km)
      61           3 :   atm::Length       topAtm(  48.0,"km");  // Upper atm. boundary for calculations
      62           3 :   atm::Pressure     Pstep(  10.0,"mb"); // Primary pressure step (10.0 mb)
      63           1 :   double   PstepFact=         1.2;     // Pressure step ratio between two consecutive layers
      64             : 
      65           2 :   atm::AtmProfile myProfile( Alt, P, T, TLR, H, WVL, Pstep, PstepFact,  topAtm, atmType );
      66             : 
      67           1 :   cout<<"BASIC ATMOSPHERIC PARAMETERS TO GENERATE REFERENCE ATMOSPHERIC PROFILE"<<endl;
      68           1 :   printf("Ground temperature T:          %.1f K \n",T.get());
      69           1 :   printf("Ground pressure P:             %.0f Pa \n",P.get("Pa")); 
      70           1 :   printf("Relative humidity rh:          %.0f percent \n",H.get("%")); 
      71           1 :   printf("Scale height h0:               %.0f m \n",WVL.get("m")); 
      72           1 :   printf("Pressure step dp:              %.0f Pa \n",Pstep.get("Pa")); 
      73           1 :   printf("Altitude alti:                 %.0f m \n",Alt.get()); 
      74           1 :   printf("Altitude top atm profile atmh: %.0f m \n",topAtm.get("m"));
      75           1 :   printf("Pressure step factordp1:       %.2f \n",PstepFact);
      76           1 :   printf("Tropospherique lapse rate:     %.1f K/km \n" , TLR);
      77             : 
      78           1 :   cout<<"Number of layers :"<<myProfile.getNumLayer()<<endl;
      79           1 :   cout<<"First guess precipitable water vapor content: " << myProfile.getGroundWH2O().get("mm") << " mm" << endl;
      80             : 
      81             :   // -------------------------------------------------------------------
      82             :   // Create a spectralGrid containing 1 astronomical band
      83             :   // This SpectralGrid is a member of a RefrativeIndexProfile object
      84             :   // -------------------------------------------------------------------
      85             : 
      86             :   // Astronomical  spectral window
      87             :   //
      88           1 :   double refFreqAstro   = 86.3*1.e9;
      89           1 :   double bandWidthAstro = 2.*1.e9; 
      90           1 :   double chanWidthAstro = 0.1*1e9; // 100Mhz
      91             :   
      92           1 :   unsigned int numChanAstro=int(bandWidthAstro/chanWidthAstro); 
      93           1 :   unsigned int refChanAstro=numChanAstro/2+1;
      94             : 
      95             :     
      96             :   atm::SpectralGrid * spGrid_p = new atm::SpectralGrid(numChanAstro,   // unsigned int numChan
      97             :                                                        refChanAstro,   // unsigned int refChan
      98           1 :                                                        atm::Frequency(refFreqAstro), // Frequency refFreq
      99           1 :                                                        atm::Frequency(chanWidthAstro));
     100             : 
     101             : 
     102             :   
     103             :     // Print description of spectral windows
     104             :     //
     105           1 :   cout<<"-----------------------------------------------"<<endl;
     106           2 :   for(unsigned int j=0; j<spGrid_p->getNumSpectralWindow(); j++){
     107           1 :     cout << "Spectral Window " << j 
     108             :       // << " Central Frequency: " <<  spGrid_p->getRefFreq(j).get("GHz") << " GHz, "
     109           2 :          << " Central Frequency: " <<  spGrid_p->getChanFreq(j,(spGrid_p->getNumChan(j)/2)).get("GHz") << " GHz, "
     110           2 :          << " Freq. Resolution: " <<  spGrid_p->getChanSep(j).get("MHz") << " MHz, "
     111           1 :          << " LO frequency: " <<  spGrid_p->getLoFrequency(j)/1.e9 << " GHz, "
     112           1 :          << " Num. of channels: " << spGrid_p->getNumChan(j)<<endl;
     113             : 
     114             :     //         for (unsigned int i=0;i<spGrid_p->getNumChan(j);i++)
     115             :     //           cout<<spGrid_p->getChanFreq(j,i).get("GHz")<<" ";
     116             :     //         cout<<endl;
     117             :   }
     118           1 :   cout<<"-----------------------------------------------"<<endl;
     119             : 
     120             : 
     121             : 
     122           1 :   cout << "Create SkyStatus object"<<endl; // and associates a WaterVaporRadiometer to it " << endl; 
     123             :         
     124             :   // Define absorption profile with SpectralGridand Profile
     125           2 :   atm::RefractiveIndexProfile absProfile(*spGrid_p, myProfile);
     126             : 
     127             :   // -------------------------------------------------------------------
     128             :   // - Create a SkyStatus object 
     129             :   // -------------------------------------------------------------------
     130             : 
     131             :   // Define skystatus instance
     132           2 :   atm::SkyStatus skyStatus(absProfile);
     133           1 :   skyStatus.setAirMass(1.);
     134           1 :   skyStatus.setUserWH2O(myProfile.getGroundWH2O());
     135             : 
     136           1 :   cout<<"end"<<endl;
     137             : 
     138             :   // -------------------------------------------------------------------
     139             :   // - Retrieve path length
     140             :   // -------------------------------------------------------------------
     141             : 
     142           1 :   double deltaPressure = 100.;
     143             :   
     144           1 :   double startPressure = 49000.;
     145           1 :   double endPressure   = 61000.;
     146           1 :   int numPressure = int((endPressure-startPressure)/deltaPressure);
     147             : 
     148           1 :   int spwid = spGrid_p->getNumSpectralWindow()-1; // last spw = astronomical band
     149             :  
     150           2 :   string filename="path.dat";
     151           1 :   ofstream resultFile(filename.c_str());
     152           1 :   resultFile <<"! Pressure TotalPath H2OPathLength DispersiveDryPathLength NonDispersiveDryPathLength"<<endl;
     153         121 :   for (unsigned int i=0; i<numPressure; i++) {
     154         120 :     skyStatus.setBasicAtmosphericParameters(atm::Pressure(startPressure+i*deltaPressure,"Pa"));
     155             :     
     156         120 :     cout<<"Number of layers :"<<skyStatus.getNumLayer()<<endl;
     157         120 :     cout<<"First guess precipitable water vapor content: " << skyStatus.getGroundWH2O().get("mm") << " mm" << endl;
     158         120 :     cout<<"Ground Pressure: " << skyStatus.getGroundPressure().get("Pa") << " Pa" << endl;
     159         120 :     cout<<"Ground Temperature: " << skyStatus.getGroundTemperature().get("K") << " K" << endl;
     160        2510 :     for (unsigned int ii=0; ii<skyStatus.getNumLayer(); ii++) {
     161        4780 :       cout << "Pressure in Layer " << ii << " = " <<  skyStatus.getLayerPressure(ii).get("Pa") << " Pa " << 
     162        2390 :         "Temperature: " <<  skyStatus.getLayerTemperature(ii).get("K") << " K " << endl;
     163             :     } 
     164             : 
     165             : 
     166         240 :     atm::Length path = skyStatus.getAverageH2OPathLength(spwid)
     167         480 :       + skyStatus.getAverageDispersiveDryPathLength(spwid)
     168         240 :       + skyStatus.getAverageNonDispersiveDryPathLength(spwid);
     169             : 
     170         120 :     resultFile<<startPressure+i*deltaPressure<<" "
     171         120 :               <<path.get()<<" "
     172         240 :               <<skyStatus.getAverageH2OPathLength(spwid).get()<<" "
     173         240 :               <<skyStatus.getAverageDispersiveDryPathLength(spwid).get()<<" "
     174         120 :               <<skyStatus.getAverageNonDispersiveDryPathLength(spwid).get()<<endl;
     175             :   }
     176           1 :   resultFile.close();
     177             : 
     178           1 :   return 0;
     179             : 
     180             : }