Main Page   Namespace List   Class Hierarchy   Compound List   File List   Namespace Members   Compound Members   File Members   Related Pages  

HeSpectrum Class Reference

Calculate the cosmic ray 4He spectrum in low earth orbit. More...

#include <HeSpectrum.h>

Inheritance diagram for HeSpectrum::

Spectrum ISpectrum List of all members.

Public Types

typedef std::pair< int, float > Intrp

Public Methods

 HeSpectrum (const std::string &params)
virtual double calculate_rate (double old_rate)
virtual double flux (double) const
 calculate flux for the current cutoff. More...

virtual double solidAngle () const
 calcualte effective solid angle for the given energy. More...

float flux (float cut) const
 Total flux in nuclei / m^2 sec ster Interpolate in table if possible, otherwise assume power law tail at high energy. More...

virtual float operator() (float) const
float cutoff () const
virtual void setPosition (double lat, double lon)
int askGPS ()
 this one asks the GPS for position. More...

float findCutoff (float rflux) const
 determine the cutoff value which will produce the desired flux. More...

float findCutoff (float lat, float lon) const
 determine the cutoff value at a geographical location. More...

virtual std::pair< float,
float > 		dir (float energy) const
 return solid angle pair (costh, phi) for the given energy. More...

virtual std::string title () const
 return a title describing the spectrum. More...

virtual const char * particleName () const
 subclasses need to specify correct particle type. More...

const char * nameOf () const
void setParticleName (std::string name)
 set the particle name. (default is "alpha"). More...


Private Methods

double HeSpectrum::flux (double lat, double lon) const
void init (float lat, float lon)
float rad2 () const
float exposure (float E) const
float cosomega (float E) const

Private Attributes

InterpVec m_en
InterpVec m_fl
InterpVec m_fluxes
float m_expo
float m_normfact
float m_tot
float m_upper
float m_etop
float m_cutoff
float m_coscutoff
float m_fluxTbl [73][13]
float m_flux
std::string m_particle_name
ObserverAdapter< HeSpectrum > m_observer

Static Private Attributes

const float m_rearth = 6371.f
const float m_altitude = 600.f

Detailed Description

Calculate the cosmic ray 4He spectrum in low earth orbit.

Uses data produced by CHIME, assuming 600 km circular orbit at solar minimum (worst case).

Interface: The constructor arguments specify the satellite position (latitude, longitude) in degrees. The position can be changed by use of the setPosition() member. The total flux in nuclei/(m^2 sec ster) is returned by the flux() member. There are 3 ways to call flux(): With no arguments it uses the current cutoff energy. If there is one argument, it is assumed to be a cutoff. Two arguments are assumed to be latitude and longitude, and the corresponding cutoff is looked up. The stored current value is not changed. Cutoff energies are returned by FindCutoff(). With no arguments, it returns the stored current value. Two arguments are assumed to be latitude and longitude, and the corresponding cutoff is looked up. The operator() function returns a sampled energy value. The argument must be a float value between 0 and 1. The dir() member returns a sampled particle's energy and direction. The argument is a value between 0 and 1, which is used to produce the energy, as in operator(). The direction is based on simple Stormer cone theory. The first component of the direction is the zenith direction cosine: +1 for particles moving downward. The second component of the direction is the earth azimuth angle: 0 from the east, +pi/2 from the north, -pi/2 from the south. At energies near the cutoff, particles come mainly from the west.
Author:
Patrick Nolan, Stanford University, 1998
$Header $

Definition at line 47 of file HeSpectrum.h.

Member Typedef Documentation

typedef std::pair<int,float> HeSpectrum::Intrp
 

Definition at line 93 of file HeSpectrum.h.

Referenced by HeSpectrum::InterpVec::interpolate(), and HeSpectrum::InterpVec::search().

Constructor & Destructor Documentation

HeSpectrum::HeSpectrum (  const std::string &    paramstring ) 
 

Definition at line 115 of file HeSpectrum.cxx.

References init(), and Spectrum::parseParamList(). 

00115                                                    {
00116     std::vector<float> params;
00117     
00118     parseParamList(paramstring,params);
00119     
00120     if(params.empty())
00121         init(0.0f,0.0f);
00122     else
00123         init(params[0], params[1]);
00124 }

Member Function Documentation

int HeSpectrum::askGPS (    ) 
 

this one asks the GPS for position.

Definition at line 190 of file HeSpectrum.cxx.

References GPS::instance(), and setPosition().

Referenced by init(). 

00191 {
00192     setPosition(GPS::instance()->lat(), GPS::instance()->lon());
00193     return 0; // can't be void in observer pattern
00194 }

double HeSpectrum::calculate_rate (  double    old_rate )  [virtual]
 

Definition at line 239 of file HeSpectrum.cxx.

References flux(), and GPS::instance(). 

00240 {
00241     return flux(GPS::instance()->lat(), GPS::instance()->lon());
00242 }

float HeSpectrum::cosomega (  float    E )  const [private]
 

Definition at line 259 of file HeSpectrum.cxx.

References m_cutoff, and rad2().

Referenced by exposure(). 

00259                                         {
00260     // Purpose:  Give opening angle of geomagnetic cutoff cone.
00261     // This is a pretty backward way to do it.  It starts from a cutoff
00262     // energy which is derived from a desired rate.  Then it calculates
00263     // the magnetic latitude from that.  Really, the cutoff and latitude
00264     // should be derived from the position by looking in tables.
00265     // Also, this is simple Størmer theory, ignoring  the penumbral
00266     // region and multipole effects.
00267     
00268     const float Mp = 4*0.938f;  // mass of alpha in GeV
00269     const float charge = 2.;    // electric charge of alpha
00270     float pcut = sqrt(m_cutoff*m_cutoff + 2.*Mp*m_cutoff);  // cutoff momentum
00271     const float moment = 59.8f; // magnetic moment of earth in convenient units
00272     float coslambda4 = 4. * pcut  * rad2() / (charge*moment);
00273     // magnetic latitude**4
00274     float p = sqrt(E*E + 2.*Mp*E); // momentum
00275     float coso = 4. * (sqrt(pcut/p) - pcut/p) * pow(coslambda4, -0.75);
00276     // opening angle of Størmer cone
00277     if (coso > 1.) return 1.;
00278     else if (coso < -1.) return -1.;
00279     else return coso;
00280 }

float HeSpectrum::cutoff (    )  const [inline]
 

Definition at line 69 of file HeSpectrum.h.

References m_cutoff. 

00069 {return m_cutoff;};

std::pair< float, float > HeSpectrum::dir (  float    energy )  const [virtual]
 

return solid angle pair (costh, phi) for the given energy.

Reimplemented from Spectrum.

Definition at line 291 of file HeSpectrum.cxx. 

00292 {
00293     
00294     // Purpose:  Random particle direction from Størmer cone
00295     
00296     // Rejection method for the direction.  Direction must be inside
00297     // the Stormer cone and above the horizon.
00298     float coszenith, sinpolar, cospolar, azi;
00299     const int try_max = 1000;
00300     static int max_tried = 0;
00301     int trial=0;
00302     do {
00303         // uniform distribution within Størmer cone
00304         cospolar = -1. + HepRandom::getTheGenerator()->flat()*
00305             (cosomega(energy)+1.);
00306         sinpolar = sqrt(1.-cospolar*cospolar);
00307         azi = 2.*M_PI* HepRandom::getTheGenerator()->flat();
00308         coszenith = cos(azi)*sinpolar;
00309     } while (coszenith < m_coscutoff && trial++ < try_max);
00310     
00311     max_tried = std::max(trial, max_tried); // keep track
00312     
00313     // Transform to local earth-based coordinates.
00314     /*float*/double earthazi = atan2(sinpolar*sin(azi), cospolar);
00315     
00316     return std::make_pair<float,float>(coszenith, earthazi);
00317 }

float HeSpectrum::exposure (  float    E )  const [private]
 

Definition at line 284 of file HeSpectrum.cxx.

References cosomega().

Referenced by setPosition(). 

00284                                         {
00285     // Geomagnetic cutoff fraction.  Varies from 0 to 1.
00286     return 0.5 * (1. + cosomega(E));
00287 }

float HeSpectrum::findCutoff (  float    lat,
float    lon
)  const
 

determine the cutoff value at a geographical location.

Definition at line 231 of file HeSpectrum.cxx.

References findCutoff(), and flux(). 

00231                                                        {
00232     // determine the cutoff value at a geographical location
00233     return findCutoff(flux(lat,lon));
00234 }

float HeSpectrum::findCutoff (  float    rflux )  const
 

determine the cutoff value which will produce the desired flux.

Definition at line 224 of file HeSpectrum.cxx.

References HeSpectrum::InterpVec::interpolate(), m_en, m_fluxes, m_upper, and HeSpectrum::InterpVec::search().

Referenced by findCutoff(), and setPosition(). 

00224                                               {
00225     // determine the cutoff value which will produce the desired flux
00226     return m_en.interpolate(m_fluxes.search(rflux-m_upper));
00227 }

float HeSpectrum::flux (  float    cut )  const
 

Total flux in nuclei / m^2 sec ster Interpolate in table if possible, otherwise assume power law tail at high energy.

Definition at line 140 of file HeSpectrum.cxx.

References HeSpectrum::InterpVec::interpolate(), m_en, m_etop, m_expo, m_fl, m_upper, and HeSpectrum::InterpVec::search(). 

00140                                       {
00141     // Purpose:  Return total flux in nuclei / m^2 sec ster
00142     // Interpolate in table if possible, otherwise assume power law
00143     //  tail at high energy.
00144     if (cut > m_etop) return m_upper * pow(cut/m_etop, m_expo);
00145     else              return m_upper + m_fl.interpolate(m_en.search(cut));
00146 }

double HeSpectrum::flux (  double    time )  const [virtual]
 

calculate flux for the current cutoff.

Reimplemented from Spectrum.

Definition at line 150 of file HeSpectrum.cxx.

References m_flux.

Referenced by calculate_rate(), findCutoff(), and setPosition(). 

00150                                     {
00151     // Purpose: calculate flux for the current cutoff
00152     return m_flux;
00153 }

double HeSpectrum::HeSpectrum::flux (  double    lat,
double    lon
)  const [private]
 

void HeSpectrum::init (  float    lat,
float    lon
)  [private]
 

Definition at line 64 of file HeSpectrum.cxx.

References askGPS(), GPS::instance(), m_altitude, m_coscutoff, m_en, m_etop, m_expo, m_fl, m_fluxes, m_fluxTbl, m_normfact, m_observer, m_particle_name, m_rearth, GPS::notification(), and setPosition().

Referenced by HeSpectrum(). 

00064                                           {
00065     //Purpose: Initializes parameters during construction
00066     int nen = sizeof(energies)/sizeof(float);
00067     m_en.reserve(nen);
00068     float amus = 4.;  // Mass of 4He; CHIME  energy is in MeV/amu
00069     int i;
00070     for (i=0; i< nen; i++) m_en.push_back(amus*energies[i]/1000.);
00071     // convert MeV to GeV
00072     
00073     m_normfact = .5*(1.+sqrt(m_altitude*(m_altitude+2.*m_rearth)) /
00074         (m_rearth+m_altitude));
00075     //geometrical shadow factor in 600 km orbit
00076     m_fluxes.reserve(nen);
00077     m_fl.reserve(nen);
00078     
00079     const float width = 0.115f; // CHIME log spacing between standard energies
00080     m_etop = m_en.back() * (1.+.5*width);  // boundary between table & tail
00081     m_expo = -2.65f + 1.0f;  // power law exponent (integral) of high-energy tail
00082     
00083     // Populate table of differential galactic fluxes -- no geomagnetic cutoff
00084     float tfl = 0.f;
00085     for (i=74; i >= 0; i--) {
00086         tfl += width*1000.*m_en[i]*m_normfact*fluxes[i];
00087         m_fluxes.insert(m_fluxes.begin(), tfl);
00088     }
00089     
00090     // table of total flux as a function of latitude and longitude
00091     for (int ii=0; ii < 73; ii++) {
00092         for (int jj=0; jj < 13; jj++) {
00093             m_fluxTbl[ii][jj] = gfluxes[jj+13*ii];
00094         }
00095     }
00096     setPosition(lat, lon);
00097     
00098     // cos(angle between zenith and horizon)
00099     m_coscutoff = -sqrt(m_altitude*m_altitude+2.*m_altitude*m_rearth)
00100         / (m_altitude+m_rearth);
00101     
00102     m_particle_name = "alpha";
00103     
00104     // set callback to be notified when the position changes
00105     m_observer.setAdapter( new ActionAdapter<HeSpectrum>(this,
00106         &HeSpectrum::askGPS) );
00107     
00108     GPS::instance()->notification().attach( &m_observer );
00109     
00110 }

const char* HeSpectrum::nameOf (    )  const [inline]
 

Definition at line 87 of file HeSpectrum.h. 

00087 {return "HeSpectrum";}

float HeSpectrum::operator() (  float    x )  const [virtual]
 

Reimplemented from Spectrum.

Definition at line 180 of file HeSpectrum.cxx. 

00180                                           {
00181     // return a random value of energy sampled from the spectrum
00182     
00183     float join = (m_tot-m_upper)/m_tot;
00184     if (x < join) return m_en.interpolate(m_fl.search((1.-x)*m_tot-m_upper));
00185     else          return m_etop*pow((1.-x)/(1.-join), 1./m_expo);
00186 }

const char * HeSpectrum::particleName (    )  const [virtual]
 

subclasses need to specify correct particle type.

Reimplemented from Spectrum.

Definition at line 130 of file HeSpectrum.cxx.

References m_particle_name. 

00130                                             {
00131     return m_particle_name.c_str();
00132 }

float HeSpectrum::rad2 (    )  const [private]
 

Definition at line 246 of file HeSpectrum.cxx.

References m_altitude, Spectrum::m_lat, Spectrum::m_lon, and m_rearth.

Referenced by cosomega(). 

00246                              {
00247     // square of (distance from center of magnetic dipole / earth radius)
00248     // Dipole is offset from the earth's center
00249     
00250     /*float*/double d2 =
00251     pow((m_rearth+m_altitude)*sin(m_lat)            - 145.1, 2) +
00252         pow((m_rearth+m_altitude)*cos(m_lat)*cos(m_lon) + 371.2, 2) +
00253         pow((m_rearth+m_altitude)*cos(m_lat)*sin(m_lon) - 233.7, 2);
00254     return d2 / (m_rearth*m_rearth);
00255 }

void HeSpectrum::setParticleName (  std::string    name ) 
 

set the particle name. (default is "alpha").

Definition at line 133 of file HeSpectrum.cxx.

References m_particle_name. 

00134 {
00135     m_particle_name = name;
00136 }

void HeSpectrum::setPosition (  double    lat,
double    lon
)  [virtual]
 

Definition at line 196 of file HeSpectrum.cxx.

References exposure(), findCutoff(), flux(), m_cutoff, m_en, m_expo, m_fl, m_flux, Spectrum::m_lat, Spectrum::m_lon, m_normfact, m_tot, and m_upper.

Referenced by askGPS(), and init(). 

00196                                                    {
00197     // Purpose:  Do the initialization necessary when moving to a new position:
00198     // look up cutoff energy, build a new table of integral alpha
00199     // fluxes
00200     
00201     m_lat = lat;
00202     m_lon = lon;
00203     
00204     // Integrated flux in the power law tail above the table.
00205     m_upper = -0.115*1000.*m_en.back()*m_normfact*fluxes[74]
00206         * pow(1.+.5*0.115,m_expo)/(0.115*m_expo);
00207     
00208     m_cutoff = findCutoff(lat,lon);
00209     
00210     // Populate table of integral fluxes modified by geomagnetic cutoff.
00211     float tfl = 0.;
00212     m_fl.erase(m_fl.begin(), m_fl.end());
00213     for (int i = 74; i >= 0; i--) {
00214         tfl += 0.115*1000.*m_en[i]*m_normfact*fluxes[i]*exposure(m_en[i]);
00215         m_fl.insert(m_fl.begin(), tfl);
00216     }
00217     
00218     m_tot = m_fl.front() + m_upper;
00219     m_flux = flux(m_cutoff);
00220 }

double HeSpectrum::solidAngle (    )  const [virtual]
 

calcualte effective solid angle for the given energy.

Reimplemented from Spectrum.

Definition at line 157 of file HeSpectrum.cxx.

References m_fluxTbl. 

00158 {
00159     // the normalization for the flux calculation
00160     return 4.*M_PI;
00161 }

std::string HeSpectrum::title (    )  const [virtual]
 

return a title describing the spectrum.

Reimplemented from Spectrum.

Definition at line 126 of file HeSpectrum.cxx. 

00126                                   {
00127     return "HeSpectrum";
00128 }

Member Data Documentation

const float HeSpectrum::m_altitude = 600.f [static, private]
 

Definition at line 61 of file HeSpectrum.cxx.

Referenced by init(), and rad2().

float HeSpectrum::m_coscutoff [private]
 

Definition at line 116 of file HeSpectrum.h.

Referenced by init().

float HeSpectrum::m_cutoff [private]
 

Definition at line 115 of file HeSpectrum.h.

Referenced by cosomega(), cutoff(), and setPosition().

InterpVec HeSpectrum::m_en [private]
 

Definition at line 107 of file HeSpectrum.h.

Referenced by findCutoff(), flux(), init(), and setPosition().

float HeSpectrum::m_etop [private]
 

Definition at line 114 of file HeSpectrum.h.

Referenced by flux(), and init().

float HeSpectrum::m_expo [private]
 

Definition at line 110 of file HeSpectrum.h.

Referenced by flux(), init(), and setPosition().

InterpVec HeSpectrum::m_fl [private]
 

Definition at line 108 of file HeSpectrum.h.

Referenced by flux(), init(), and setPosition().

float HeSpectrum::m_flux [private]
 

Definition at line 118 of file HeSpectrum.h.

Referenced by flux(), and setPosition().

InterpVec HeSpectrum::m_fluxes [private]
 

Definition at line 109 of file HeSpectrum.h.

Referenced by findCutoff(), and init().

float HeSpectrum::m_fluxTbl[73][13] [private]
 

Definition at line 117 of file HeSpectrum.h.

Referenced by init(), and solidAngle().

float HeSpectrum::m_normfact [private]
 

Definition at line 111 of file HeSpectrum.h.

Referenced by init(), and setPosition().

ObserverAdapter< HeSpectrum > HeSpectrum::m_observer [private]
 

Definition at line 127 of file HeSpectrum.h.

Referenced by init().

std::string HeSpectrum::m_particle_name [private]
 

Definition at line 126 of file HeSpectrum.h.

Referenced by init(), particleName(), and setParticleName().

const float HeSpectrum::m_rearth = 6371.f [static, private]
 

Definition at line 60 of file HeSpectrum.cxx.

Referenced by init(), and rad2().

float HeSpectrum::m_tot [private]
 

Definition at line 112 of file HeSpectrum.h.

Referenced by setPosition().

float HeSpectrum::m_upper [private]
 

Definition at line 113 of file HeSpectrum.h.

Referenced by findCutoff(), flux(), and setPosition().

The documentation for this class was generated from the following files:
Generated on Wed Oct 16 14:01:34 2002 by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001