physics_funcs.cc

Go to the documentation of this file.

/* Copyright (C) 2002-2008
   Patrick Eriksson <Patrick.Eriksson@rss.chalmers.se>
   Stefan Buehler   <sbuehler@ltu.se>
 
   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version.
 
   This program 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 General Public License for more details.
 
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
   USA. */
 
 
 
 
/*===========================================================================
  === File description 
  ===========================================================================*/
 
/*===========================================================================
  === External declarations
  ===========================================================================*/
 
#include <cmath>
#include <stdexcept>
#include "physics_funcs.h"
#include "messages.h"          
#include "mystring.h"
#include "physics_funcs.h"
 
extern const Numeric BOLTZMAN_CONST;
extern const Numeric DEG2RAD;
extern const Numeric PLANCK_CONST;
extern const Numeric SPEED_OF_LIGHT;
 
 
 
/*===========================================================================
  === The functions (in alphabetical order)
  ===========================================================================*/
 
 
Numeric dinvplanckdI(
        const Numeric&  i,
        const Numeric&  f )
{
  assert( i >= 0 );
  assert( f > 0 );
 
  // Use always double to avoid numerical problem (see invrayjean)
  static const double a = PLANCK_CONST / BOLTZMAN_CONST;
  static const double b = 2 * PLANCK_CONST / ( SPEED_OF_LIGHT*SPEED_OF_LIGHT );
 
  const double d    = b * pow((double)f,3) / (double)i;
  const double binv = a * (double)f / log( d + 1 );
 
  return pow(binv,2) / ( a * (double)f * (double)i * (1/d+1) ); 
}
 
 
 
 
void fresnel(
             Complex&   Rv,
             Complex&   Rh,
       const Complex&   n1,
       const Complex&   n2,
       const Numeric&   theta )
{
  const Numeric theta1 = DEG2RAD * theta;
  const Numeric costheta1 = cos( theta1 );
  const Numeric costheta2 = cos( asin( n1.real() * sin(theta1) / n2.real() ) );
 
  Complex a, b;
  a  = n2 * costheta1;
  b  = n1 * costheta2;
  Rv = ( a - b ) / ( a + b );
  a  = n1 * costheta1;
  b  = n2 * costheta2;
  Rh = ( a - b ) / ( a + b );
}
 
 
 
Numeric invplanck(
        const Numeric&  i,
        const Numeric&  f )
{
  assert( i >= 0 );
  assert( f > 0 );
 
  // Use always double to avoid numerical problem (see invrayjean)
  static const double a = PLANCK_CONST / BOLTZMAN_CONST;
  // Double clealrly needed. b likely too small for float. b is approx 1e-50
  static const double b = 2 * PLANCK_CONST / ( SPEED_OF_LIGHT*SPEED_OF_LIGHT );
 
  return   ( a * (double)f ) / log( (b*pow((double)f,3))/(double)i + 1 );
}
 
 
 
 
Numeric invrayjean(
        const Numeric&  i,
        const Numeric&  f )
{
  assert( f > 0 );
 
  // Double must be used here (if not, the result can be NaN when using float)
  
  static const double   a = SPEED_OF_LIGHT*SPEED_OF_LIGHT/(2*BOLTZMAN_CONST);
 
  return  ( a * (double)i ) / ( (double)f * (double)f );
}
 
 
 
Numeric number_density(  
        const Numeric&   p,
        const Numeric&   t )
{
  assert( p >= 0 );
  assert( t > 0 );
  return   p / ( t * BOLTZMAN_CONST );
}
 
 
 
 
Numeric planck( 
        const Numeric&   f, 
        const Numeric&   t )
{
  assert( t >= 0 );
  assert( f > 0 );
 
  if( t == 0 )
    { return 0; }
 
  else
    {
    // Double must be used here (if not, a becomes 0 when using float)
    static const double a = 2 * PLANCK_CONST / (SPEED_OF_LIGHT*SPEED_OF_LIGHT);
    static const double b = PLANCK_CONST / BOLTZMAN_CONST;
  
    return  ( a * pow((double)f,3) ) / ( exp( (b*(double)f)/(double)t ) - 1 );
    }
}
 
 
 
 
Numeric rayjean(
        const Numeric&  f,
        const Numeric&  t )
{
  assert( f > 0 );
 
  // Double must be used here (if not, the result can be NaN when using float)
  
  static const double   a = SPEED_OF_LIGHT*SPEED_OF_LIGHT/(2*BOLTZMAN_CONST);
 
  return  ( (double)f * (double)f ) / ( a * (double)t );
}