casa6

/*******************************************************************************

 * ALMA - Atacama Large Millimiter Array

 * (c) Instituto de Estructura de la Materia, 2009

 *

 * This library is free software; you can redistribute it and/or

 * modify it under the terms of the GNU Lesser General Public

 * License as published by the Free Software Foundation; either

 * version 2.1 of the License, or (at your option) any later version.

 *

 * This library 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

 * Lesser General Public License for more details.

 *

 * You should have received a copy of the GNU Lesser General Public

 * License along with this library; if not, write to the Free Software

 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA

 *

 * "@(#) $Id: ATMRefractiveIndex.cpp Exp $"

 *

 * who       when      what

 * --------  --------  ----------------------------------------------

 * pardo     24/03/09  created

 */

​

#include "ATMRefractiveIndex.h"

​

#include <stdio.h>

#include <iostream>

#include <math.h>

​

using namespace std;

​

using std::cout;

​

​

ATM_NAMESPACE_BEGIN

​

  // Constructors

​

​

  RefractiveIndex::RefractiveIndex(){ }

​

  RefractiveIndex::~RefractiveIndex()

  {

    // rmRefractiveIndex();  // nonexistent method

  }

​

  std::complex<double> RefractiveIndex::getRefractivity_o2(double temperature, double pressure, double wvpressure, double frequency){

​

    static const double abun_18o=0.0020439;

    static const double abun_17o=0.0003750;

    static const double o2_mixing_ratio=0.2092;

​

    std::complex<double> ccc = (mkSpecificRefractivity_16o16o(temperature,pressure,wvpressure,frequency)*(1.0-2.0*(abun_18o+abun_17o))*

               (1.0-exp(-1556.38*1.43/temperature))+

               mkSpecificRefractivity_16o16o_vib(temperature,pressure,wvpressure,frequency)*(1.0-2.0*(abun_18o+abun_17o))*

               exp(-1556.38*1.43/temperature)

               +mkSpecificRefractivity_16o18o(temperature,pressure,wvpressure,frequency)*2.0*abun_18o

               +mkSpecificRefractivity_16o17o(temperature,pressure,wvpressure,frequency)*2.0*abun_17o

               )*o2_mixing_ratio*pressure*100.0/(1.380662e-23*temperature);

​

    //      if(frequency<143&&frequency>142.21){cout << "O2: " << frequency << "  " << ccc << "  " << pressure <<  endl;}

​

    return ccc;

​

  }

​

  std::complex<double> RefractiveIndex::getRefractivity_o2(double temperature,double pressure,double wvpressure,

                              double frequency,double width,unsigned int n)

  {

    std::complex<double> average(0.0,0.0);

    double newfreq;

    for(unsigned int i=0; i<n; i++){

      if(n==1){

    newfreq=frequency;

      }else{

    newfreq=frequency-(width/2.0)+(width/(n-1))*i;

      }

      average=average+getRefractivity_o2(temperature,pressure,wvpressure,newfreq);

    }

    std::complex<double> averagen(real(average)/n,imag(average)/n);

    return averagen;

  }

​

​

​

  std::complex<double> RefractiveIndex::getRefractivity_h2o(double temperature, double pressure, double wvpressure, double frequency){

​

    static const double abun_18o=0.0020439;

    static const double abun_17o=0.0003750;