Line data    Source code

       1             : /** 
       2             :    Bojan Nikolic <b.nikolic@mrao.cam.ac.uk>, <bojan@bnikolic.co.uk>
       3             :    Initial version June 2009
       4             :    Maintained by ESO since 2013.
       5             : 
       6             :    This file is part of LibAIR and is licensed under GNU Public
       7             :    License Version 2
       8             : 
       9             :    \file model_water.hpp
      10             :    Renamed model_water.h 2023
      11             : 
      12             :  */
      13             : #ifndef _LIBAIR_MODEL_WATER_HPP__
      14             : #define _LIBAIR_MODEL_WATER_HPP__
      15             : 
      16             : #include <memory>
      17             : 
      18             : #include "model_iface.h"
      19             : #include "radiometermeasure.h"
      20             : #include "basicphys.h"
      21             : #include "cloudywater.h"
      22             : 
      23             : namespace LibAIR2 {
      24             : 
      25             :   // Forward decleration
      26             :   class Radiometer;
      27             :   
      28             :   template <class AM>
      29             :   class WaterModel:
      30             :     public WVRAtmoQuantModel
      31             :   {
      32             : 
      33             :     std::shared_ptr<Radiometer> mr;
      34             :     std::shared_ptr<AM> am;
      35             : 
      36             :   public:
      37             : 
      38             :     // -------------------- Public data ---------------------
      39             : 
      40             :     inline double n_bump( ) const;
      41             :     inline double tau_bump( ) const;
      42             : 
      43             :     // -------------------- Construction / Destruction ------
      44             : 
      45             : 
      46          81 :     WaterModel(std::shared_ptr<Radiometer> mr,
      47             :                std::shared_ptr<AM> am):
      48          81 :       mr(mr),
      49          81 :       am(am)
      50             :     {
      51          81 :     }
      52             :     
      53             :     // Inherited from WVRAtmoQuants
      54           0 :     double eval(size_t ch) const
      55             :     {
      56           0 :       return mr->eval(am->TbGrid(), 
      57           0 :                       ch);
      58             :     }
      59             : 
      60    79505305 :     virtual void eval(std::vector<double> & res) const
      61             :     {
      62    79505305 :       const size_t nc = mr->nchannels();
      63    79505305 :       res.resize(nc);
      64             :       const std::vector<double> & 
      65    79505305 :         tbg = am->TbGrid();
      66             :       
      67   397526525 :       for (size_t i = 0 ; i < nc ; ++i)
      68   318021220 :         res[i]= mr->eval(tbg,i);
      69    79505305 :     }
      70             : 
      71           0 :     virtual double dTdc (size_t ch) const
      72             :     {
      73           0 :       const double invdelta = 0.5/n_bump( );
      74           0 :       const double on=am->n;
      75             :       
      76           0 :       am->n=on+n_bump( );
      77           0 :       const double fwdval =  eval(ch);
      78           0 :       am->n=on-n_bump( );
      79           0 :       const double backval = eval(ch);
      80           0 :       am->n=on;
      81             :       
      82           0 :       return (fwdval - backval) * invdelta ;
      83             :     }
      84             : 
      85      807800 :     void dTdc (std::vector<double> &res) const
      86             :     {
      87      807800 :       const double invdelta = 0.5/ n_bump( );
      88      807800 :       const size_t nc = mr->nchannels();
      89      807800 :       const double on=am->n;
      90             :       
      91      807800 :       std::vector<double> fwdval(nc);
      92      807800 :       am->n=on+n_bump( );
      93      807800 :       eval(fwdval);
      94             :       
      95             : 
      96      807800 :       std::vector<double> backval(nc);
      97      807800 :       am->n=on-n_bump( );
      98      807800 :       eval(backval);
      99             : 
     100      807800 :       am->n=on;
     101             :       
     102      807800 :       res.resize(nc);
     103     4039000 :       for(size_t  i=0 ; i < res.size() ; ++i)
     104     3231200 :         res[i] = (fwdval[i] - backval[i]) * invdelta ;
     105             :       
     106      807800 :     }
     107             : 
     108           0 :     double dTdL_ND (size_t ch) const
     109             :     {
     110           0 :       return dTdc(ch) / SW_WaterToPath_Simplified(1.0, 
     111           0 :                                                   am->T);
     112             :     }
     113             : 
     114      807800 :     void dTdL_ND (std::vector<double> & res) const
     115             :     {
     116      807800 :       dTdc(res);
     117      807800 :       const double conv = 1.0 / SW_WaterToPath_Simplified(1.0, 
     118      807800 :                                                           am->T);
     119             :       
     120     4039000 :       for (size_t i =0 ; i < res.size() ; ++i)
     121     3231200 :         res[i] *= conv;
     122      807800 :     }
     123             : 
     124             :     /** Compute the derivative of sky temperatures wvr the continuum
     125             :         opacity term
     126             :      */
     127             :     void dTdTau (std::vector<double> &) const
     128             :     {
     129             :       throw std::runtime_error("This function requires specific impelmentation classes");
     130             :     }    
     131             : 
     132             :     // Inherited from WVRAtmoModel
     133         243 :     void AddParams ( std::vector< Minim::DParamCtr > &pars )
     134             :     {
     135         243 :       pars.push_back(Minim::DParamCtr ( &am->n ,      
     136             :                                         "n", 
     137             :                                         true     ,                       
     138             :                                         "Water column (mm)"
     139             :                                         ));
     140             :       
     141         243 :       pars.push_back(Minim::DParamCtr ( &am->T ,      
     142             :                                         "T", 
     143             :                                         true     ,                       
     144             :                                         "Temperature (K)"
     145             :                                         ));
     146             :       
     147         243 :       pars.push_back(Minim::DParamCtr ( &am->P ,      
     148             :                                         "P", 
     149             :                                         true     ,                       
     150             :                                         "Pressure (mBar)"
     151             :                                         ));    
     152             : 
     153         243 :     }
     154             : 
     155             :   };
     156             : 
     157             :   template<>   inline 
     158           0 :   void WaterModel<ICloudyWater>::AddParams (std::vector< Minim::DParamCtr > &pars)
     159             :   {
     160           0 :       pars.push_back(Minim::DParamCtr ( &am->n ,      
     161             :                                         "n", 
     162             :                                         true     ,                       
     163             :                                         "Water column (mm)"
     164             :                                         ));
     165             :       
     166           0 :       pars.push_back(Minim::DParamCtr ( &am->T ,      
     167             :                                         "T", 
     168             :                                         true     ,                       
     169             :                                         "Temperature (K)"
     170             :                                         ));
     171             :       
     172           0 :       pars.push_back(Minim::DParamCtr ( &am->P ,      
     173             :                                         "P", 
     174             :                                         true     ,                       
     175             :                                         "Pressure (mBar)"
     176             :                                         ));    
     177             : 
     178           0 :       pars.push_back(Minim::DParamCtr ( &am->tau183 ,      
     179             :                                         "tau183", 
     180             :                                         true     ,                       
     181             :                                         "Opacity at 183GHz"
     182             :                                         ));    
     183             : 
     184           0 :   }
     185             : 
     186             :   template<>   inline double WaterModel<ISingleLayerWater>::tau_bump( ) const { return 0.001; }
     187     2423400 :   template<>   inline double WaterModel<ISingleLayerWater>::n_bump( ) const { return 0.001; }
     188             :   template<>   inline double WaterModel<ICloudyWater>::tau_bump( ) const { return 0.001; }
     189           0 :   template<>   inline double WaterModel<ICloudyWater>::n_bump( ) const { return 0.001; }
     190             : 
     191             :   template<>   inline 
     192             :   void WaterModel<ICloudyWater>::dTdTau (std::vector<double> &res) const
     193             :   {
     194             :     const double invdelta = 0.5/ tau_bump( );
     195             :     const size_t nc = mr->nchannels();
     196             :     const double on=am->tau183;
     197             :     
     198             :     std::vector<double> fwdval(nc);
     199             :     am->tau183=on+tau_bump( );
     200             :     eval(fwdval);
     201             :     
     202             :     
     203             :     std::vector<double> backval(nc);
     204             :     am->tau183=on-tau_bump( );
     205             :     eval(backval);
     206             :     
     207             :     am->tau183=on;
     208             :     
     209             :     res.resize(nc);
     210             :     for(size_t  i=0 ; i < res.size() ; ++i)
     211             :       res[i] = (fwdval[i] - backval[i]) * invdelta ;
     212             :     
     213             :   }    
     214             : 
     215             : }
     216             : 
     217             : 
     218             : 
     219             : #endif