text:
#ifndef __INC_ESTIMATE_THREADED__
#define __INC_ESTIMATE_THREADED__

#include <unistd.h>
#include <pthread.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

/* common data block for most options to be passed around */

struct optstruct;
struct modelstruct;

//! used to pass the args into estimate_thread_function
struct estimate_thetas_params{
	// these are direct copies of the respective data structures
	struct optstruct* options;
	struct modelstruct* the_model;
	// things needed specifically for the estimation
	gsl_rng* random_number;
	gsl_matrix* h_matrix;
	int max_tries;
	int success_count;
	double lhood_current;
	double my_best; // want to check on the local best values
} estimate_thetas_params;

#define USEMUTEX

#include "../modelstruct.h"
#include "../optstruct.h"

void estimate_thetas_threaded(modelstruct* the_model, optstruct* options);

// this won't work unless it has access to the nasty globals in estimate_threaded.c
void* estimate_thread_function(void* args);
int get_number_cpus(void);

void setup_params(struct estimate_thetas_params *params, modelstruct* the_model, optstruct* options, int nthreads, int max_tries);

void fprintPt(FILE *f, pthread_t pt);

// see the source for defs of the number of threads etc etc
#endif

max_stars_repo_name: jackdawjackdaw/emulator

max_issues_repo_name: MADAI/MADAIEmulator

max_forks_repo_name: MADAI/MADAIEmulator


text:
/* Author:  G. Jungman */

/* Convenience header */
#ifndef __GSL_SPECFUNC_H__
#define __GSL_SPECFUNC_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_sf.h>

#endif /* __GSL_SPECFUNC_H__ */

ext: h

max_stars_repo_name: berkus/music-cs

max_issues_repo_name: berkus/music-cs

max_forks_repo_name: berkus/music-cs


text:
#include <stdio.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_randist.h>

int
dofit(const gsl_multifit_robust_type *T,
      const gsl_matrix *X, const gsl_vector *y,
      gsl_vector *c, gsl_matrix *cov)
{
  int s;
  gsl_multifit_robust_workspace * work
    = gsl_multifit_robust_alloc (T, X->size1, X->size2);

  s = gsl_multifit_robust (X, y, c, cov, work);
  gsl_multifit_robust_free (work);

  return s;
}

int
main (int argc, char **argv)
{
  int i;
  size_t n;
  const size_t p = 2; /* linear fit */
  gsl_matrix *X, *cov;
  gsl_vector *x, *y, *c, *c_ols;
  const double a = 1.45; /* slope */
  const double b = 3.88; /* intercept */
  gsl_rng *r;

  if (argc != 2)
    {
      fprintf (stderr,"usage: robfit n\n");
      exit (-1);
    }

  n = atoi (argv[1]);

  X = gsl_matrix_alloc (n, p);
  x = gsl_vector_alloc (n);
  y = gsl_vector_alloc (n);

  c = gsl_vector_alloc (p);
  c_ols = gsl_vector_alloc (p);
  cov = gsl_matrix_alloc (p, p);

  r = gsl_rng_alloc(gsl_rng_default);

  /* generate linear dataset */
  for (i = 0; i < n - 3; i++)
    {
      double dx = 10.0 / (n - 1.0);
      double ei = gsl_rng_uniform(r);
      double xi = -5.0 + i * dx;
      double yi = a * xi + b;

      gsl_vector_set (x, i, xi);
      gsl_vector_set (y, i, yi + ei);
    }

  /* add a few outliers */
  gsl_vector_set(x, n - 3, 4.7);
  gsl_vector_set(y, n - 3, -8.3);

  gsl_vector_set(x, n - 2, 3.5);
  gsl_vector_set(y, n - 2, -6.7);

  gsl_vector_set(x, n - 1, 4.1);
  gsl_vector_set(y, n - 1, -6.0);

  /* construct design matrix X for linear fit */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);

      gsl_matrix_set (X, i, 0, 1.0);
      gsl_matrix_set (X, i, 1, xi);
    }

  /* perform robust and OLS fit */
  dofit(gsl_multifit_robust_ols, X, y, c_ols, cov);
  dofit(gsl_multifit_robust_bisquare, X, y, c, cov);

  /* output data and model */
  for (i = 0; i < n; ++i)
    {
      double xi = gsl_vector_get(x, i);
      double yi = gsl_vector_get(y, i);
      gsl_vector_view v = gsl_matrix_row(X, i);
      double y_ols, y_rob, y_err;

      gsl_multifit_robust_est(&v.vector, c, cov, &y_rob, &y_err);
      gsl_multifit_robust_est(&v.vector, c_ols, cov, &y_ols, &y_err);

      printf("%g %g %g %g\n", xi, yi, y_rob, y_ols);
    }

#define C(i) (gsl_vector_get(c,(i)))
#define COV(i,j) (gsl_matrix_get(cov,(i),(j)))

  {
    printf ("# best fit: Y = %g + %g X\n",
            C(0), C(1));

    printf ("# covariance matrix:\n");
    printf ("# [ %+.5e, %+.5e\n",
               COV(0,0), COV(0,1));
    printf ("#   %+.5e, %+.5e\n",
               COV(1,0), COV(1,1));
  }

  gsl_matrix_free (X);
  gsl_vector_free (x);
  gsl_vector_free (y);
  gsl_vector_free (c);
  gsl_vector_free (c_ols);
  gsl_matrix_free (cov);
  gsl_rng_free(r);

  return 0;
}

ext: c

max_stars_repo_name: ruslankuzmin/julia

max_issues_repo_name: ruslankuzmin/julia

max_forks_repo_name: ruslankuzmin/julia


text:
/* sum/demo.c
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * 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 3 of the License, 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#include <stdio.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sum.h>

#define N 20

int
main (void)
{
  double t[N];
  double sum_accel, err;
  double sum = 0;
  int n;

  gsl_sum_levin_u_workspace * w = gsl_sum_levin_u_alloc (N);

  const double zeta_2 = M_PI * M_PI / 6.0;

  /* terms for zeta(2) = \sum_{n=1}^{\infty} 1/n^2 */

  for (n = 0; n < N; n++)
    {
      double np1 = n + 1.0;
      t[n] = 1.0 / (np1 * np1);
      sum += t[n] ;
    }

  gsl_sum_levin_u_accel (t, N, w, &sum_accel, &err);

  printf("term-by-term sum = % .16f using %d terms\n", sum, N) ;

  printf("term-by-term sum = % .16f using %d terms\n",
         w->sum_plain, w->terms_used) ;

  printf("exact value      = % .16f\n", zeta_2) ;
  printf("accelerated sum  = % .16f using %d terms\n",
         sum_accel, w->terms_used) ;

  printf("estimated error  = % .16f\n", err) ;
  printf("actual error     = % .16f\n", sum_accel - zeta_2) ;

  gsl_sum_levin_u_free (w);

  return 0;
}

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
// Copyright (c) Microsoft Corporation.
// Licensed under the MIT License.

#pragma once

#include <memory>
#include <future>
#include <array>
#include <gsl\gsl>

#include <boost\optional.hpp>

#include "Data\Data.h"
#include "MageSettings.h"
#include "Device\Device.h"

#include <SensorSample.h>
#include <Plat/CameraDevice/CameraSettings.h>

namespace mage
{
    class PoseHistory;

    /* Entry point of the MAGESlam library */
    class MAGESlam
    {
        struct Impl;
    public:

        struct CameraConfiguration
        {
            CameraIdentity CameraIdentity;
            Size Size{};
            PixelFormat Format;
            Matrix Extrinsics;

            CameraConfiguration(mage::CameraIdentity cameraIdentity, const mage::Size& size, PixelFormat format, const Matrix& extrinsics)
                : CameraIdentity{ cameraIdentity },
                Format{ format },
                Size{ size },
                Extrinsics{ extrinsics }
            {
            };
        };

        struct FrameFormat
        {
            const mage::FrameId FrameId;
            std::shared_ptr<const calibration::CameraModel> CameraModel;
            std::chrono::system_clock::time_point Timestamp;
            mira::CameraSettings CameraSettings;

            FrameFormat(
                const mage::FrameId& frameId,
                std::shared_ptr<const calibration::CameraModel> cameraModel,
                const std::chrono::system_clock::time_point& timestamp,
                const mira::CameraSettings& cameraSettings)
                : FrameId{ frameId },
                CameraModel{ cameraModel },
                Timestamp{ timestamp },
                CameraSettings{ cameraSettings }
            {}
        };

        struct Tracking
        {
            /*
                The Matrix is a row-major column-vector view matrix. right handed
                    R, R, R, Tx
                    R, R, R, Ty
                    R, R, R, Tz
                    0, 0, 0, 1
            */
            Matrix Pose;
            TrackingState State;

            Tracking(const Matrix& mat, TrackingState state)
                : Pose{ mat }, State{ state }
            {}

            Tracking()
                : Pose{ }, State{ TrackingState::SKIPPED }
            {}

            bool IsPoseGood() const
            {
                return State == TrackingState::TRACKING;
            }

            bool IsPoseSkipped() const
            {
                return State == TrackingState::SKIPPED;
            }
        };

        struct TrackedFrame
        {
            Tracking Tracking;
            std::shared_ptr<const calibration::CameraModel> CameraModel;
            Depth Depth;

            TrackedFrame(const MAGESlam::Tracking& tracking, std::shared_ptr<const calibration::CameraModel> cameraModel, const mage::Depth& depth)
                : Tracking{ tracking }, CameraModel{ cameraModel }, Depth{ depth }
            {}

            TrackedFrame() = default;
        };

        class FossilizedMap
        {
        public:
            ~FossilizedMap();

            /*
                Get the most up-to-date pose for the specified frame ID
            */
            std::vector<boost::optional<MAGESlam::TrackedFrame>> GetTrackingResultsForFrames(const gsl::span<const FrameId> frameIds) const;

            bool TryGetVolumeOfInterest(AxisAlignedVolume& volumeOfInterest) const;

        private:
            FossilizedMap(std::unique_ptr<const PoseHistory> history, std::vector<Position> mapPoints, const MageSlamSettings& settings);

            friend class MAGESlam;

            struct Impl;
            std::unique_ptr<Impl> m_impl;
        };

        explicit MAGESlam(const MageSlamSettings&, gsl::span<const CameraConfiguration> cameras, const device::IMUCharacterization& imu);

        MAGESlam(const MAGESlam&) = delete;
        MAGESlam& operator =(const MAGESlam&) = delete;

        struct Frame
        {
            FrameFormat Format;
            gsl::span<const uint8_t> Bytes;

            Frame(const FrameFormat& format, gsl::span<const uint8_t> bytes)
                : Format(format),
                Bytes(bytes)
            {}
        };

        /*
            Processes an image frame and returns a future Pose and the state of the
            tracker while it was generating it.
        */
        std::future<Tracking> ProcessFrame(const Frame& frame);

        /*
            Processes an stereo image pair and returns a future Pose for each frame and the state of the
            tracker while it was generating it.
        */
        std::pair<std::future<Tracking>, std::future<Tracking>> ProcessStereoFrames(const Frame& one, const Frame& two);

        /*
            Get the most up-to-date pose for the specified frame ID
        */
        std::vector<boost::optional<MAGESlam::TrackedFrame>> GetTrackingResultsForFrames(const gsl::span<const FrameId> frameIds) const;

        // adds a sensor sample to the sample queue for sensor fusion
        void AddSensorSample(const SensorSample& sample);

        // get normalized vector in mage's world coordinate frame describing the fuser's gravity estimate
        // returns whether the gravity estimate is valid.
        bool GetGravityDirection(Direction& gravDir) const;

        // gets the scale to apply to mage to bring it to meters as estimated by
        // the stereo tracking code
        float GetStereoMageMeterEstimate() const;

        // gets the scale to apply to mage to bring it to meters as estimated by the fuser
        // returns whether the scale is valid
        bool  GetScaleFromIMU(float& scaleMAGEtoMeters) const;

        bool TryGetVolumeOfInterest(AxisAlignedVolume& volumeOfInterest) const;

        // Destroys the current instance of MAGESlam and converts it to a fossilized map
        // that can be used for offline queries of frame poses
        static std::unique_ptr<FossilizedMap> Fossilize(std::unique_ptr<MAGESlam> slam);

        ~MAGESlam();
    private:
        std::unique_ptr<Impl> m_impl;
    };
}

max_stars_repo_name: syntheticmagus/mageslam

max_issues_repo_name: syntheticmagus/mageslam

max_forks_repo_name: syntheticmagus/mageslam


text:
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_poly.h>
#include "userFunc.h"

#define G        6.67384e-8    /* 2010 CODATA value in CGS units */
#define C        2.99792458e10 /* Exact value in CGS units */
#define ALPHAMIN 1.0e-3        /* Minimum alpha allowed for numerical reasons */
#define YR       (365.25*24.0*3600.0)  /* 1 yr in sec */
#define MYR      (1.0e6*YR)    /* 1 Myr in seconds */
#define MSUN     1.9886e33     /* Solar mass in g */

/**************************************************************************/
/* This defines userFunc routines for the Krumholz & Kruijssen CMZ        */
/* disk model including feedback                                          */
/**************************************************************************/

/* Define the parameters to be used here */
typedef struct {
  double eta;
  double sigmath;
  double shapefac;
  double zetad;
  double alphacoef;
  double epsff0;
  long tctr;
  double p_per_sn;
  double e_per_sn;
  double star_rad_var;
  double col_min;
  double sigma_max;
  double col_floor;
  double en_wind;
  void *wksp;   /* This a placeholder for the extra memory stored at
		   the end of the parameters for temporary workspace
		   and to store the SF history. The size of storage
		   here will be (nbin_SB99+5)*grd->nr; we implement
		   pointers to this data via macro, below */
} paramType;

#define MSRC_TMP(p)     ((double *) &(((paramType *) (p))->wksp))
#define PDOT(p,g)       ((double *) &(((paramType *) (p))->wksp) + (g)->nr)
#define EDOT(p,g)       ((double *) &(((paramType *) (p))->wksp) + 2*(g)->nr)
#define MDOT(p,g)       ((double *) &(((paramType *) (p))->wksp) + 3*(g)->nr)
#define SIGMASTAR(p,g)  ((double *) &(((paramType *) (p))->wksp) + 4*(g)->nr)

/* Declare some utility functions to compute the SFR, mass loss rate,
   and momentum and energy injection rate per unit area */
void injection_rate(double *mstar, double p_per_sn, double e_per_sn,
		    double *pdot_tot, double *edot_tot, double *mdot_tot);
void sfrwind(const grid *grd, const double *col, const double *pres,
	     void *params, double *sigma_SFR, double *sigma_wind);

/* Data from starburst99 */
const int nbin_SB99  = 1000;      /* Number of output times */
const double dt_SB99 = 0.05;      /* Output time interval in Myr */

/* Supernova rate, expressed as log SNe / year / 10^6 Msun of stars, for stellar populations of different ages */
const double snr_SB99[1000] =
  {  -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00,  -3.64,  -3.33,  -3.32,  -3.31,  -3.31,  -3.28,  -3.30,  -3.27,  -3.29,  -3.28,  -3.28,  -3.25,  -3.25,  -3.27,  -3.26,  -3.24,  -3.27,  -3.24,  -3.26,  -3.24,  -3.26,  -3.26,  -3.26,  -3.25,  -3.25,  -3.25,  -3.25,  -3.25,  -3.25,  -3.22,  -3.22,  -3.24,  -3.24,  -3.22,  -3.22,  -3.24,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.24,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.26,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.27,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.28,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.29,  -3.30,  -3.30,  -3.30,  -3.33,  -3.33,  -3.33,  -3.34,  -3.34,  -3.34,  -3.34,  -3.35,  -3.35,  -3.35,  -3.35,  -3.35,  -3.35,  -3.36,  -3.36,  -3.36,  -3.36,  -3.36,  -3.37,  -3.37,  -3.37,  -3.37,  -3.37,  -3.37,  -3.38,  -3.38,  -3.38,  -3.38,  -3.38,  -3.38,  -3.39,  -3.39,  -3.39,  -3.39,  -3.39,  -3.39,  -3.40,  -3.40,  -3.40,  -3.40,  -3.40,  -3.40,  -3.40,  -3.41,  -3.41,  -3.41,  -3.41,  -3.41,  -3.41,  -3.42,  -3.42,  -3.42,  -3.42,  -3.42,  -3.42,  -3.42,  -3.42,  -3.43,  -3.43,  -3.43,  -3.43,  -3.43,  -3.43,  -3.43,  -3.44,  -3.44,  -3.44,  -3.44,  -3.44,  -3.44,  -3.44,  -3.44,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.45,  -3.46,  -3.46,  -3.46,  -3.46,  -3.46,  -3.46,  -3.46,  -3.46,  -3.47,  -3.47,  -3.47,  -3.47,  -3.47,  -3.47,  -3.47,  -3.47,  -3.47,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.48,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.49,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.50,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.51,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.52,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.53,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.54,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.55,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.56,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.57,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.58,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.59,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.60,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.61,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.62,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.63,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.64,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.65,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.66,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.70,  -3.68,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.67,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.68,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.69,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.70,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.71,  -3.72,  -3.72,  -3.72,  -3.72,  -3.72,  -3.72,  -3.72,  -3.74, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00, -30.00 };

/* Luminosity of stellar populations of different ages, expressed as log of the luminosity in erg/s per 10^6 Msun worth of stars */
const double lum_SB99[1000] =
  {  42.496, 42.497, 42.498, 42.500, 42.501, 42.503, 42.505, 42.506, 42.508, 42.510, 42.512, 42.513, 42.515, 42.517, 42.519, 42.520, 42.522, 42.524, 42.526, 42.528, 42.530, 42.531, 42.533, 42.535, 42.537, 42.539, 42.541, 42.543, 42.545, 42.547, 42.549, 42.551, 42.553, 42.556, 42.558, 42.560, 42.562, 42.564, 42.567, 42.569, 42.572, 42.574, 42.576, 42.579, 42.582, 42.584, 42.587, 42.590, 42.593, 42.595, 42.598, 42.601, 42.603, 42.606, 42.609, 42.611, 42.614, 42.618, 42.619, 42.619, 42.616, 42.612, 42.607, 42.602, 42.593, 42.579, 42.566, 42.554, 42.541, 42.529, 42.517, 42.505, 42.493, 42.482, 42.471, 42.460, 42.450, 42.440, 42.431, 42.422, 42.413, 42.405, 42.396, 42.388, 42.380, 42.372, 42.364, 42.357, 42.349, 42.341, 42.334, 42.326, 42.319, 42.312, 42.304, 42.297, 42.289, 42.282, 42.274, 42.266, 42.258, 42.250, 42.243, 42.236, 42.228, 42.221, 42.214, 42.207, 42.200, 42.193, 42.187, 42.180, 42.174, 42.167, 42.161, 42.155, 42.149, 42.143, 42.136, 42.130, 42.124, 42.118, 42.113, 42.107, 42.101, 42.095, 42.089, 42.084, 42.078, 42.072, 42.067, 42.061, 42.055, 42.050, 42.044, 42.038, 42.033, 42.027, 42.021, 42.015, 42.010, 42.004, 41.999, 41.994, 41.988, 41.983, 41.978, 41.973, 41.967, 41.962, 41.957, 41.952, 41.947, 41.943, 41.937, 41.932, 41.927, 41.922, 41.917, 41.913, 41.908, 41.903, 41.898, 41.894, 41.889, 41.884, 41.879, 41.874, 41.870, 41.865, 41.860, 41.855, 41.851, 41.846, 41.841, 41.837, 41.833, 41.829, 41.824, 41.821, 41.817, 41.813, 41.809, 41.805, 41.801, 41.797, 41.794, 41.790, 41.786, 41.783, 41.779, 41.776, 41.772, 41.768, 41.765, 41.762, 41.758, 41.755, 41.751, 41.748, 41.745, 41.742, 41.738, 41.735, 41.732, 41.729, 41.726, 41.723, 41.719, 41.716, 41.713, 41.710, 41.707, 41.704, 41.701, 41.698, 41.695, 41.693, 41.690, 41.687, 41.684, 41.681, 41.679, 41.676, 41.673, 41.671, 41.667, 41.665, 41.662, 41.659, 41.656, 41.654, 41.651, 41.648, 41.645, 41.642, 41.639, 41.636, 41.634, 41.631, 41.628, 41.625, 41.622, 41.620, 41.617, 41.614, 41.611, 41.608, 41.606, 41.603, 41.599, 41.597, 41.594, 41.591, 41.588, 41.586, 41.583, 41.580, 41.578, 41.575, 41.572, 41.569, 41.566, 41.564, 41.561, 41.559, 41.556, 41.553, 41.551, 41.548, 41.545, 41.543, 41.540, 41.538, 41.535, 41.533, 41.530, 41.528, 41.525, 41.523, 41.521, 41.518, 41.516, 41.513, 41.511, 41.509, 41.506, 41.504, 41.502, 41.499, 41.497, 41.495, 41.492, 41.490, 41.488, 41.486, 41.484, 41.482, 41.480, 41.478, 41.475, 41.473, 41.471, 41.469, 41.467, 41.464, 41.462, 41.460, 41.458, 41.456, 41.454, 41.452, 41.450, 41.448, 41.445, 41.444, 41.442, 41.440, 41.438, 41.436, 41.434, 41.432, 41.430, 41.428, 41.426, 41.425, 41.423, 41.421, 41.419, 41.417, 41.416, 41.414, 41.412, 41.410, 41.409, 41.407, 41.406, 41.404, 41.402, 41.401, 41.399, 41.397, 41.396, 41.394, 41.392, 41.391, 41.390, 41.388, 41.387, 41.385, 41.384, 41.382, 41.381, 41.379, 41.378, 41.376, 41.375, 41.374, 41.373, 41.371, 41.370, 41.368, 41.367, 41.365, 41.364, 41.363, 41.361, 41.360, 41.359, 41.357, 41.356, 41.355, 41.353, 41.352, 41.351, 41.349, 41.348, 41.347, 41.346, 41.344, 41.343, 41.342, 41.340, 41.339, 41.338, 41.336, 41.335, 41.334, 41.332, 41.331, 41.330, 41.329, 41.328, 41.326, 41.325, 41.324, 41.323, 41.322, 41.320, 41.319, 41.318, 41.317, 41.315, 41.314, 41.313, 41.312, 41.311, 41.309, 41.308, 41.307, 41.306, 41.305, 41.304, 41.303, 41.302, 41.301, 41.300, 41.298, 41.297, 41.296, 41.295, 41.294, 41.292, 41.291, 41.290, 41.289, 41.288, 41.287, 41.286, 41.285, 41.284, 41.283, 41.282, 41.281, 41.279, 41.279, 41.277, 41.277, 41.276, 41.274, 41.273, 41.272, 41.271, 41.270, 41.269, 41.268, 41.267, 41.266, 41.265, 41.264, 41.263, 41.262, 41.261, 41.259, 41.259, 41.258, 41.257, 41.256, 41.254, 41.253, 41.252, 41.252, 41.251, 41.249, 41.248, 41.247, 41.246, 41.245, 41.244, 41.243, 41.242, 41.241, 41.240, 41.239, 41.238, 41.238, 41.236, 41.235, 41.234, 41.233, 41.232, 41.232, 41.231, 41.230, 41.230, 41.229, 41.227, 41.227, 41.225, 41.225, 41.224, 41.223, 41.222, 41.221, 41.220, 41.219, 41.218, 41.218, 41.217, 41.216, 41.215, 41.214, 41.213, 41.212, 41.211, 41.211, 41.209, 41.208, 41.208, 41.207, 41.206, 41.205, 41.204, 41.203, 41.203, 41.202, 41.201, 41.200, 41.199, 41.198, 41.198, 41.197, 41.196, 41.195, 41.194, 41.193, 41.192, 41.191, 41.191, 41.190, 41.189, 41.188, 41.187, 41.187, 41.186, 41.185, 41.184, 41.184, 41.183, 41.182, 41.181, 41.180, 41.180, 41.179, 41.178, 41.177, 41.176, 41.175, 41.175, 41.174, 41.173, 41.172, 41.171, 41.171, 41.170, 41.169, 41.168, 41.167, 41.166, 41.165, 41.165, 41.164, 41.163, 41.162, 41.162, 41.161, 41.160, 41.159, 41.159, 41.158, 41.157, 41.156, 41.156, 41.155, 41.154, 41.153, 41.152, 41.151, 41.150, 41.150, 41.149, 41.148, 41.147, 41.146, 41.146, 41.145, 41.144, 41.143, 41.143, 41.141, 41.141, 41.140, 41.139, 41.139, 41.138, 41.137, 41.136, 41.136, 41.135, 41.134, 41.134, 41.133, 41.118, 41.131, 41.130, 41.130, 41.129, 41.128, 41.128, 41.127, 41.126, 41.125, 41.125, 41.124, 41.123, 41.122, 41.121, 41.120, 41.119, 41.119, 41.118, 41.117, 41.116, 41.115, 41.114, 41.113, 41.113, 41.112, 41.111, 41.111, 41.110, 41.109, 41.108, 41.107, 41.107, 41.106, 41.105, 41.104, 41.104, 41.103, 41.102, 41.101, 41.100, 41.099, 41.099, 41.098, 41.097, 41.096, 41.095, 41.095, 41.094, 41.093, 41.093, 41.092, 41.091, 41.090, 41.089, 41.088, 41.088, 41.087, 41.086, 41.085, 41.085, 41.084, 41.084, 41.083, 41.082, 41.082, 41.081, 41.080, 41.079, 41.079, 41.078, 41.077, 41.076, 41.076, 41.075, 41.074, 41.073, 41.073, 41.072, 41.071, 41.070, 41.070, 41.069, 41.068, 41.068, 41.067, 41.066, 41.065, 41.065, 41.064, 41.063, 41.063, 41.062, 41.061, 41.060, 41.060, 41.059, 41.058, 41.058, 41.057, 41.056, 41.055, 41.055, 41.054, 41.053, 41.053, 41.052, 41.051, 41.050, 41.050, 41.049, 41.048, 41.048, 41.047, 41.046, 41.045, 41.045, 41.044, 41.043, 41.043, 41.042, 41.041, 41.041, 41.040, 41.040, 41.039, 41.038, 41.037, 41.036, 41.036, 41.035, 41.035, 41.034, 41.033, 41.033, 41.032, 41.031, 41.030, 41.030, 41.029, 41.028, 41.028, 41.027, 41.027, 41.026, 41.025, 41.024, 41.024, 41.023, 41.023, 41.022, 41.021, 41.020, 41.020, 41.019, 41.018, 41.017, 41.017, 41.016, 41.016, 41.015, 41.014, 41.014, 41.013, 41.012, 41.012, 41.011, 41.010, 41.010, 41.009, 41.008, 41.008, 41.007, 41.007, 41.006, 41.005, 41.005, 41.004, 41.003, 41.003, 41.002, 41.001, 41.000, 41.000, 40.999, 40.999, 40.998, 40.998, 40.997, 40.997, 40.996, 40.996, 40.995, 40.994, 40.994, 40.993, 40.993, 40.992, 40.992, 40.991, 40.990, 40.990, 40.989, 40.988, 40.988, 40.987, 40.986, 40.985, 40.985, 40.984, 40.984, 40.983, 40.983, 40.982, 40.981, 40.981, 40.980, 40.979, 40.979, 40.978, 40.978, 40.977, 40.976, 40.976, 40.975, 40.975, 40.974, 40.973, 40.973, 40.972, 40.971, 40.971, 40.970, 40.970, 40.969, 40.969, 40.968, 40.967, 40.967, 40.966, 40.966, 40.965, 40.964, 40.964, 40.963, 40.962, 40.962, 40.962, 40.961, 40.960, 40.942, 40.959, 40.958, 40.958, 40.957, 40.957, 40.956, 40.955, 40.955, 40.954, 40.954, 40.953, 40.953, 40.952, 40.951, 40.951, 40.950, 40.949, 40.949, 40.948, 40.948, 40.947, 40.946, 40.946, 40.945, 40.945, 40.944, 40.943, 40.943, 40.942, 40.942, 40.941, 40.941, 40.940, 40.939, 40.939, 40.938, 40.938, 40.937, 40.937, 40.936, 40.935, 40.935, 40.934, 40.934, 40.933, 40.933, 40.932, 40.931, 40.931, 40.930, 40.929, 40.929, 40.928, 40.928, 40.928, 40.927, 40.926, 40.926, 40.925, 40.925, 40.924, 40.923, 40.923, 40.922, 40.922, 40.921, 40.921, 40.920, 40.920, 40.919, 40.919, 40.918, 40.918, 40.917, 40.916, 40.916, 40.915, 40.915, 40.914, 40.914, 40.913, 40.912, 40.912, 40.911, 40.911, 40.911, 40.910, 40.909, 40.909, 40.908, 40.907, 40.907, 40.906, 40.906, 40.905, 40.905, 40.904, 40.904, 40.903, 40.902, 40.902, 40.901, 40.901, 40.900, 40.899, 40.899, 40.899, 40.898, 40.898, 40.897, 40.897, 40.897, 40.896, 40.895, 40.895, 40.894, 40.894, 40.893, 40.893, 40.892, 40.892, 40.891, 40.891, 40.890, 40.890, 40.889, 40.889, 40.888, 40.888, 40.887, 40.887, 40.886, 40.886, 40.885, 40.884, 40.884, 40.883, 40.883, 40.882, 40.881, 40.881, 40.880, 40.880, 40.879, 40.879, 40.878, 40.877 };

/* Wind momentum flux, expressed as log of the momentum in dyne per 10^6 Msun worth of stars */
const double pdot_wind_SB99[1000] =
  {  31.689, 31.695, 31.700, 31.705, 31.709, 31.714, 31.718, 31.723, 31.728, 31.732, 31.737, 31.741, 31.746, 31.750, 31.755, 31.759, 31.764, 31.768, 31.773, 31.778, 31.782, 31.787, 31.792, 31.796, 31.801, 31.805, 31.809, 31.814, 31.818, 31.822, 31.826, 31.829, 31.833, 31.836, 31.838, 31.841, 31.842, 31.843, 31.844, 31.843, 31.842, 31.839, 31.835, 31.830, 31.822, 31.812, 31.803, 31.821, 31.844, 31.868, 31.893, 31.929, 31.982, 32.001, 32.016, 32.012, 32.033, 32.064, 32.168, 32.222, 32.272, 32.308, 32.325, 32.339, 32.340, 32.336, 32.332, 32.329, 32.329, 32.295, 32.251, 32.214, 32.183, 32.160, 32.140, 32.118, 32.088, 32.011, 31.980, 31.946, 31.916, 31.885, 31.868, 31.842, 31.808, 31.805, 31.747, 31.751, 31.724, 31.707, 31.702, 31.673, 31.674, 31.630, 31.612, 31.595, 31.552, 31.545, 31.543, 31.526, 31.508, 31.494, 31.476, 31.421, 31.395, 31.389, 31.375, 31.373, 31.385, 31.157, 31.352, 31.337, 31.317, 31.352, 31.273, 31.245, 31.226, 31.185, 31.173, 31.125, 31.050, 31.063, 31.102, 31.113, 31.140, 31.160, 31.155, 31.046, 31.019, 31.004, 30.984, 30.981, 30.978, 30.979, 30.964, 30.954, 30.947, 30.937, 30.894, 30.895, 30.893, 30.896, 30.895, 30.895, 30.897, 30.898, 30.896, 30.895, 30.900, 30.903, 30.905, 30.905, 30.907, 30.909, 30.910, 30.914, 30.916, 30.911, 30.907, 30.892, 30.878, 30.859, 30.839, 30.818, 30.799, 30.777, 30.759, 30.739, 30.720, 30.699, 30.681, 30.664, 30.642, 30.624, 30.608, 30.585, 30.570, 30.549, 30.533, 30.514, 30.499, 30.478, 30.461, 30.444, 30.431, 30.414, 30.396, 30.383, 30.362, 30.351, 30.336, 30.315, 30.304, 30.285, 30.277, 30.260, 30.242, 30.235, 30.213, 30.203, 30.191, 30.179, 30.160, 30.149, 30.138, 30.124, 30.113, 30.094, 30.082, 30.073, 30.063, 30.043, 30.033, 30.023, 30.012, 29.999, 29.983, 29.972, 29.962, 29.951, 29.941, 29.930, 29.913, 29.903, 29.891, 29.878, 29.871, 29.861, 29.854, 29.846, 29.841, 29.821, 29.813, 29.798, 29.786, 29.775, 29.758, 29.745, 29.722, 29.700, 29.689, 29.650, 29.644, 29.604, 29.584, 29.571, 29.557, 29.544, 29.532, 29.520, 29.509, 29.500, 29.491, 29.482, 29.474, 29.464, 29.456, 29.447, 29.437, 29.429, 29.422, 29.415, 29.408, 29.401, 29.394, 29.387, 29.380, 29.373, 29.366, 29.360, 29.353, 29.346, 29.340, 29.333, 29.326, 29.320, 29.313, 29.307, 29.300, 29.293, 29.287, 29.281, 29.274, 29.268, 29.261, 29.255, 29.249, 29.242, 29.236, 29.230, 29.224, 29.217, 29.211, 29.205, 29.198, 29.192, 29.186, 29.180, 29.174, 29.168, 29.162, 29.156, 29.150, 29.144, 29.138, 29.131, 29.125, 29.119, 29.113, 29.107, 29.101, 29.096, 29.090, 29.084, 29.078, 29.072, 29.066, 29.061, 29.055, 29.049, 29.044, 29.038, 29.032, 29.026, 29.021, 29.015, 29.009, 29.004, 28.999, 28.993, 28.988, 28.982, 28.977, 28.972, 28.966, 28.961, 28.956, 28.951, 28.945, 28.941, 28.936, 28.930, 28.925, 28.920, 28.915, 28.910, 28.905, 28.900, 28.895, 28.889, 28.884, 28.879, 28.874, 28.868, 28.863, 28.858, 28.853, 28.848, 28.842, 28.837, 28.832, 28.827, 28.822, 28.817, 28.812, 28.807, 28.802, 28.797, 28.792, 28.787, 28.782, 28.777, 28.772, 28.768, 28.764, 28.758, 28.753, 28.749, 28.744, 28.739, 28.735, 28.730, 28.726, 28.721, 28.716, 28.712, 28.708, 28.703, 28.699, 28.694, 28.690, 28.685, 28.680, 28.677, 28.672, 28.667, 28.663, 28.659, 28.655, 28.651, 28.647, 28.642, 28.638, 28.634, 28.630, 28.626, 28.622, 28.618, 28.614, 28.609, 28.606, 28.602, 28.598, 28.594, 28.590, 28.587, 28.582, 28.578, 28.575, 28.571, 28.567, 28.563, 28.560, 28.556, 28.552, 28.548, 28.545, 28.541, 28.537, 28.534, 28.530, 28.527, 28.524, 28.519, 28.516, 28.513, 28.509, 28.506, 28.502, 28.499, 28.495, 28.491, 28.488, 28.484, 28.481, 28.478, 28.475, 28.472, 28.469, 28.465, 28.461, 28.458, 28.455, 28.447, 28.448, 28.445, 28.442, 28.439, 28.435, 28.433, 28.429, 28.426, 28.422, 28.419, 28.416, 28.413, 28.407, 28.406, 28.404, 28.400, 28.397, 28.395, 28.392, 28.386, 28.383, 28.385, 28.377, 28.374, 28.374, 28.369, 28.365, 28.366, 28.362, 28.360, 28.354, 28.355, 28.350, 28.347, 28.341, 28.342, 28.339, 28.339, 28.334, 28.331, 28.328, 28.326, 28.324, 28.321, 28.318, 28.315, 28.312, 28.310, 28.307, 28.305, 28.410, 28.373, 28.447, 28.424, 28.464, 28.443, 28.488, 28.484, 28.489, 28.487, 28.488, 28.486, 28.492, 28.488, 28.493, 28.491, 28.488, 28.485, 28.482, 28.486, 28.480, 28.482, 28.483, 28.478, 28.476, 28.478, 28.476, 28.475, 28.476, 28.469, 28.471, 28.466, 28.466, 28.470, 28.468, 28.462, 28.460, 28.462, 28.459, 28.456, 28.453, 28.454, 28.455, 28.453, 28.449, 28.448, 28.444, 28.444, 28.440, 28.441, 28.437, 28.435, 28.433, 28.436, 28.434, 28.431, 28.429, 28.426, 28.425, 28.423, 28.421, 28.417, 28.416, 28.415, 28.412, 28.409, 28.408, 28.405, 28.405, 28.401, 28.400, 28.397, 28.394, 28.394, 28.390, 28.388, 28.384, 28.383, 28.380, 28.376, 28.375, 28.372, 28.371, 28.371, 28.368, 28.366, 28.363, 28.363, 28.361, 28.359, 28.358, 28.358, 28.353, 28.354, 28.352, 28.351, 28.350, 28.350, 28.136, 28.349, 28.348, 28.346, 28.346, 28.346, 28.345, 28.343, 28.340, 28.340, 28.340, 28.339, 28.339, 28.336, 28.334, 28.333, 28.332, 28.331, 28.329, 28.329, 28.328, 28.328, 28.325, 28.323, 28.324, 28.321, 28.319, 28.318, 28.318, 28.317, 28.316, 28.315, 28.314, 28.312, 28.311, 28.310, 28.310, 28.307, 28.306, 28.305, 28.304, 28.303, 28.301, 28.301, 28.300, 28.299, 28.298, 28.297, 28.298, 28.294, 28.293, 28.294, 28.293, 28.291, 28.291, 28.289, 28.289, 28.288, 28.286, 28.285, 28.283, 28.284, 28.283, 28.283, 28.282, 28.280, 28.279, 28.277, 28.278, 28.276, 28.274, 28.273, 28.273, 28.273, 28.273, 28.270, 28.269, 28.269, 28.267, 28.266, 28.265, 28.266, 28.265, 28.263, 28.261, 28.261, 28.261, 28.259, 28.259, 28.258, 28.256, 28.256, 28.254, 28.254, 28.251, 28.253, 28.252, 28.252, 28.250, 28.249, 28.246, 28.248, 28.245, 28.245, 28.244, 28.244, 28.240, 28.239, 28.239, 28.238, 28.239, 28.238, 28.237, 28.234, 28.235, 28.234, 28.235, 28.235, 28.233, 28.229, 28.230, 28.226, 28.229, 28.229, 28.226, 28.225, 28.224, 28.225, 28.224, 28.223, 28.223, 28.220, 28.219, 28.218, 28.220, 28.218, 28.217, 28.216, 28.213, 28.213, 28.214, 28.212, 28.214, 28.209, 28.214, 28.208, 28.210, 28.205, 28.205, 28.206, 28.208, 28.203, 28.203, 28.202, 28.202, 28.201, 28.201, 28.201, 28.200, 28.198, 28.194, 28.193, 28.193, 28.199, 28.196, 28.196, 28.195, 28.193, 28.195, 28.189, 28.189, 28.188, 28.188, 28.186, 28.189, 28.186, 28.188, 28.187, 28.188, 28.181, 28.181, 28.181, 28.178, 28.178, 28.178, 28.176, 28.176, 28.176, 28.174, 28.175, 28.174, 28.173, 28.173, 28.172, 28.170, 28.170, 28.172, 28.171, 28.168, 28.168, 28.169, 28.167, 28.168, 28.167, 28.167, 28.164, 28.166, 28.162, 28.165, 28.165, 28.158, 28.157, 28.157, 28.155, 28.159, 28.156, 28.155, 28.156, 28.160, 28.151, 28.155, 28.149, 28.149, 28.155, 28.148, 28.148, 28.152, 28.145, 28.150, 28.150, 28.146, 28.146, 28.144, 28.143, 28.143, 28.142, 28.143, 28.142, 28.142, 28.142, 28.142, 28.136, 28.141, 28.139, 28.138, 28.135, 28.135, 28.138, 28.135, 28.129, 27.728, 28.133, 28.132, 28.132, 28.129, 28.131, 28.131, 28.129, 28.123, 28.122, 28.124, 28.121, 28.119, 28.123, 28.118, 28.125, 28.125, 28.116, 28.119, 28.119, 28.113, 28.114, 28.114, 28.114, 28.115, 28.112, 28.113, 28.110, 28.113, 28.108, 28.114, 28.106, 28.106, 28.112, 28.110, 28.112, 28.111, 28.111, 28.109, 28.102, 28.102, 28.098, 28.108, 28.109, 28.100, 28.107, 28.102, 28.106, 28.105, 28.105, 28.096, 28.095, 28.091, 28.095, 28.092, 28.092, 28.092, 28.098, 28.089, 28.091, 28.098, 28.098, 28.090, 28.090, 28.086, 28.086, 28.096, 28.087, 28.092, 28.092, 28.092, 28.092, 28.088, 28.087, 28.084, 28.082, 28.084, 28.088, 28.088, 28.088, 28.088, 28.088, 28.076, 28.087, 28.087, 28.079, 28.087, 28.073, 28.078, 28.083, 28.083, 28.073, 28.081, 28.072, 28.080, 28.071, 28.073, 28.073, 28.072, 28.076, 28.080, 28.072, 28.066, 28.073, 28.078, 28.078, 28.071, 28.071, 28.071, 28.066, 28.066, 28.065, 28.056, 28.069, 28.053, 28.051, 28.061, 28.049, 28.049, 28.049, 28.051, 28.048, 28.048, 28.041, 28.041, 28.041, 28.040, 28.040, 28.041, 28.039, 28.039, 28.026, 28.035, 28.026, 28.026, 28.025, 28.025, 28.032, 28.023, 28.031, 28.021, 28.021, 28.021, 28.020, 28.017, 28.015, 28.017, 28.016 };

/* Wind energy flux, expressed as log of the power in erg/s per 10^6 Msun worth of stars */
const double edot_wind_SB99[1000] =
  { 39.937, 39.941, 39.944, 39.948, 39.951, 39.953, 39.957, 39.960, 39.963, 39.966, 39.969, 39.972, 39.975, 39.977, 39.980, 39.983, 39.986, 39.989, 39.991, 39.994, 39.997, 39.999, 40.002, 40.004, 40.006, 40.008, 40.010, 40.012, 40.013, 40.014, 40.015, 40.016, 40.016, 40.016, 40.015, 40.014, 40.012, 40.009, 40.006, 40.001, 39.995, 39.988, 39.980, 39.970, 39.958, 39.943, 39.929, 39.937, 39.949, 39.960, 39.971, 39.991, 40.022, 40.026, 40.026, 40.012, 40.012, 40.011, 40.116, 40.169, 40.215, 40.250, 40.273, 40.298, 40.303, 40.295, 40.288, 40.278, 40.272, 40.243, 40.208, 40.179, 40.154, 40.132, 40.113, 40.087, 40.055, 39.949, 39.898, 39.849, 39.803, 39.750, 39.722, 39.691, 39.649, 39.644, 39.587, 39.595, 39.564, 39.545, 39.541, 39.510, 39.512, 39.453, 39.429, 39.418, 39.393, 39.385, 39.379, 39.350, 39.323, 39.296, 39.260, 39.211, 39.181, 39.169, 39.145, 39.138, 39.146, 38.963, 39.105, 39.091, 39.070, 39.126, 39.017, 38.986, 38.962, 38.918, 38.901, 38.856, 38.788, 38.790, 38.812, 38.819, 38.833, 38.843, 38.837, 38.762, 38.743, 38.723, 38.701, 38.690, 38.684, 38.683, 38.670, 38.663, 38.656, 38.650, 38.629, 38.629, 38.627, 38.629, 38.629, 38.630, 38.633, 38.635, 38.635, 38.635, 38.639, 38.642, 38.645, 38.646, 38.649, 38.650, 38.655, 38.659, 38.662, 38.659, 38.653, 38.638, 38.621, 38.602, 38.580, 38.557, 38.537, 38.514, 38.494, 38.473, 38.453, 38.432, 38.413, 38.396, 38.374, 38.356, 38.341, 38.320, 38.306, 38.285, 38.271, 38.253, 38.239, 38.219, 38.204, 38.188, 38.175, 38.160, 38.142, 38.130, 38.111, 38.101, 38.087, 38.068, 38.058, 38.041, 38.033, 38.017, 38.001, 37.994, 37.975, 37.965, 37.954, 37.942, 37.926, 37.915, 37.905, 37.892, 37.883, 37.865, 37.854, 37.845, 37.835, 37.818, 37.808, 37.798, 37.789, 37.776, 37.762, 37.751, 37.741, 37.731, 37.721, 37.711, 37.695, 37.685, 37.673, 37.660, 37.653, 37.643, 37.635, 37.627, 37.620, 37.601, 37.593, 37.578, 37.566, 37.555, 37.537, 37.525, 37.503, 37.482, 37.471, 37.436, 37.429, 37.393, 37.373, 37.360, 37.347, 37.334, 37.322, 37.310, 37.299, 37.289, 37.279, 37.270, 37.261, 37.251, 37.242, 37.233, 37.223, 37.214, 37.206, 37.198, 37.191, 37.183, 37.175, 37.167, 37.160, 37.152, 37.145, 37.137, 37.129, 37.122, 37.114, 37.107, 37.099, 37.091, 37.084, 37.077, 37.069, 37.061, 37.054, 37.047, 37.039, 37.032, 37.024, 37.017, 37.009, 37.002, 36.995, 36.987, 36.980, 36.973, 36.965, 36.957, 36.950, 36.943, 36.936, 36.928, 36.921, 36.914, 36.906, 36.899, 36.891, 36.884, 36.877, 36.869, 36.862, 36.855, 36.847, 36.839, 36.832, 36.825, 36.817, 36.810, 36.803, 36.796, 36.789, 36.782, 36.775, 36.769, 36.762, 36.755, 36.748, 36.741, 36.735, 36.728, 36.720, 36.715, 36.708, 36.701, 36.695, 36.688, 36.682, 36.675, 36.668, 36.663, 36.656, 36.650, 36.643, 36.637, 36.631, 36.624, 36.618, 36.612, 36.606, 36.599, 36.593, 36.587, 36.581, 36.575, 36.569, 36.563, 36.557, 36.550, 36.545, 36.539, 36.533, 36.527, 36.520, 36.515, 36.509, 36.503, 36.497, 36.491, 36.486, 36.480, 36.474, 36.468, 36.463, 36.456, 36.451, 36.445, 36.440, 36.434, 36.429, 36.423, 36.417, 36.412, 36.406, 36.400, 36.395, 36.390, 36.383, 36.378, 36.373, 36.368, 36.362, 36.357, 36.352, 36.346, 36.341, 36.336, 36.330, 36.325, 36.320, 36.314, 36.309, 36.304, 36.299, 36.294, 36.288, 36.283, 36.278, 36.273, 36.268, 36.263, 36.258, 36.253, 36.248, 36.242, 36.238, 36.233, 36.228, 36.223, 36.219, 36.213, 36.208, 36.203, 36.199, 36.194, 36.189, 36.184, 36.180, 36.175, 36.170, 36.165, 36.160, 36.155, 36.151, 36.146, 36.142, 36.136, 36.133, 36.126, 36.122, 36.118, 36.113, 36.109, 36.104, 36.099, 36.095, 36.090, 36.085, 36.080, 36.076, 36.072, 36.067, 36.063, 36.059, 36.054, 36.049, 36.045, 36.039, 36.036, 36.030, 36.027, 36.022, 36.019, 36.013, 36.009, 36.005, 36.000, 35.996, 35.991, 35.987, 35.982, 35.979, 35.973, 35.969, 35.965, 35.960, 35.957, 35.952, 35.948, 35.944, 35.940, 35.936, 35.930, 35.926, 35.922, 35.917, 35.914, 35.909, 35.905, 35.901, 35.897, 35.893, 35.890, 35.884, 35.879, 35.875, 35.870, 35.867, 35.863, 35.859, 35.855, 35.851, 35.846, 35.842, 35.838, 35.834, 35.830, 35.824, 35.822, 35.978, 35.924, 36.029, 35.997, 36.053, 36.025, 36.084, 36.079, 36.087, 36.083, 36.086, 36.083, 36.091, 36.086, 36.094, 36.091, 36.087, 36.084, 36.080, 36.087, 36.079, 36.082, 36.084, 36.076, 36.075, 36.078, 36.076, 36.075, 36.076, 36.067, 36.069, 36.062, 36.062, 36.068, 36.064, 36.057, 36.053, 36.056, 36.052, 36.048, 36.043, 36.045, 36.047, 36.043, 36.039, 36.038, 36.031, 36.032, 36.028, 36.029, 36.023, 36.020, 36.019, 36.022, 36.020, 36.016, 36.015, 36.010, 36.009, 36.007, 36.004, 35.999, 35.997, 35.997, 35.992, 35.989, 35.988, 35.985, 35.985, 35.980, 35.978, 35.975, 35.971, 35.970, 35.965, 35.962, 35.958, 35.957, 35.953, 35.948, 35.947, 35.944, 35.943, 35.942, 35.940, 35.937, 35.933, 35.934, 35.932, 35.930, 35.928, 35.929, 35.923, 35.924, 35.922, 35.921, 35.921, 35.921, 35.395, 35.919, 35.917, 35.916, 35.914, 35.913, 35.912, 35.911, 35.907, 35.907, 35.906, 35.905, 35.906, 35.902, 35.900, 35.898, 35.897, 35.896, 35.893, 35.894, 35.892, 35.892, 35.888, 35.886, 35.887, 35.885, 35.880, 35.880, 35.879, 35.879, 35.878, 35.877, 35.874, 35.873, 35.872, 35.870, 35.871, 35.868, 35.866, 35.866, 35.864, 35.862, 35.860, 35.860, 35.859, 35.858, 35.858, 35.856, 35.857, 35.852, 35.852, 35.852, 35.852, 35.848, 35.850, 35.847, 35.847, 35.846, 35.842, 35.842, 35.839, 35.841, 35.839, 35.840, 35.839, 35.837, 35.836, 35.834, 35.835, 35.831, 35.831, 35.829, 35.830, 35.830, 35.829, 35.825, 35.825, 35.824, 35.822, 35.821, 35.821, 35.822, 35.821, 35.819, 35.817, 35.816, 35.815, 35.814, 35.813, 35.813, 35.810, 35.811, 35.808, 35.808, 35.805, 35.807, 35.806, 35.806, 35.803, 35.803, 35.800, 35.802, 35.799, 35.798, 35.798, 35.798, 35.792, 35.791, 35.791, 35.791, 35.792, 35.790, 35.790, 35.786, 35.788, 35.787, 35.789, 35.788, 35.786, 35.781, 35.781, 35.777, 35.781, 35.781, 35.777, 35.777, 35.775, 35.777, 35.776, 35.774, 35.775, 35.772, 35.770, 35.769, 35.772, 35.769, 35.769, 35.768, 35.764, 35.763, 35.766, 35.763, 35.766, 35.759, 35.766, 35.758, 35.762, 35.755, 35.755, 35.757, 35.760, 35.753, 35.753, 35.753, 35.753, 35.752, 35.752, 35.752, 35.752, 35.748, 35.743, 35.742, 35.742, 35.751, 35.747, 35.747, 35.747, 35.744, 35.747, 35.739, 35.739, 35.738, 35.738, 35.736, 35.740, 35.735, 35.739, 35.739, 35.740, 35.730, 35.731, 35.731, 35.727, 35.727, 35.727, 35.725, 35.725, 35.725, 35.723, 35.724, 35.723, 35.723, 35.723, 35.721, 35.719, 35.719, 35.722, 35.721, 35.717, 35.717, 35.719, 35.717, 35.719, 35.717, 35.717, 35.714, 35.717, 35.712, 35.717, 35.716, 35.705, 35.705, 35.705, 35.703, 35.710, 35.705, 35.705, 35.706, 35.712, 35.700, 35.706, 35.697, 35.697, 35.706, 35.697, 35.697, 35.703, 35.694, 35.701, 35.700, 35.696, 35.696, 35.693, 35.693, 35.693, 35.692, 35.693, 35.692, 35.692, 35.692, 35.692, 35.684, 35.692, 35.689, 35.688, 35.683, 35.683, 35.688, 35.685, 35.677, 34.086, 35.682, 35.683, 35.682, 35.679, 35.682, 35.682, 35.679, 35.671, 35.671, 35.674, 35.670, 35.668, 35.675, 35.667, 35.678, 35.679, 35.666, 35.670, 35.671, 35.663, 35.665, 35.666, 35.666, 35.667, 35.663, 35.666, 35.662, 35.667, 35.660, 35.669, 35.658, 35.659, 35.668, 35.664, 35.667, 35.668, 35.668, 35.665, 35.656, 35.656, 35.651, 35.665, 35.667, 35.655, 35.666, 35.658, 35.664, 35.664, 35.664, 35.650, 35.651, 35.644, 35.651, 35.647, 35.647, 35.647, 35.657, 35.644, 35.647, 35.657, 35.658, 35.647, 35.647, 35.642, 35.642, 35.657, 35.645, 35.652, 35.652, 35.653, 35.653, 35.648, 35.646, 35.642, 35.640, 35.643, 35.649, 35.649, 35.649, 35.650, 35.650, 35.634, 35.650, 35.650, 35.639, 35.651, 35.631, 35.638, 35.645, 35.645, 35.632, 35.644, 35.631, 35.643, 35.631, 35.634, 35.634, 35.633, 35.639, 35.644, 35.634, 35.626, 35.636, 35.644, 35.644, 35.635, 35.635, 35.635, 35.627, 35.628, 35.628, 35.614, 35.633, 35.610, 35.607, 35.623, 35.606, 35.606, 35.606, 35.609, 35.606, 35.606, 35.596, 35.596, 35.596, 35.596, 35.595, 35.597, 35.595, 35.594, 35.577, 35.589, 35.577, 35.577, 35.576, 35.576, 35.587, 35.575, 35.586, 35.573, 35.573, 35.573, 35.572, 35.570, 35.567, 35.571, 35.570 };

/* Cumulative mass lost, expressed as log of the mass in Msun per 10^6
   Msun of stars */
const double mass_return_SB99[1000] =
  {  1.743, 2.225, 2.451, 2.600, 2.713, 2.803, 2.879, 2.944, 3.001, 3.053, 3.099, 3.142, 3.181, 3.218, 3.252, 3.284, 3.314, 3.343, 3.370, 3.396, 3.421, 3.445, 3.468, 3.491, 3.512, 3.533, 3.553, 3.572, 3.591, 3.610, 3.628, 3.645, 3.662, 3.679, 3.695, 3.711, 3.727, 3.742, 3.756, 3.771, 3.785, 3.799, 3.812, 3.825, 3.837, 3.849, 3.860, 3.872, 3.885, 3.898, 3.913, 3.928, 3.946, 3.965, 3.985, 4.004, 4.024, 4.053, 4.088, 4.123, 4.159, 4.196, 4.230, 4.261, 4.296, 4.333, 4.366, 4.396, 4.425, 4.450, 4.472, 4.491, 4.508, 4.523, 4.537, 4.551, 4.563, 4.575, 4.586, 4.596, 4.606, 4.616, 4.626, 4.635, 4.645, 4.654, 4.663, 4.672, 4.682, 4.691, 4.700, 4.709, 4.718, 4.726, 4.734, 4.742, 4.749, 4.756, 4.763, 4.770, 4.777, 4.783, 4.789, 4.795, 4.800, 4.806, 4.812, 4.817, 4.822, 4.828, 4.833, 4.837, 4.842, 4.846, 4.850, 4.854, 4.858, 4.862, 4.866, 4.870, 4.873, 4.877, 4.880, 4.884, 4.887, 4.891, 4.894, 4.897, 4.900, 4.904, 4.907, 4.910, 4.913, 4.917, 4.920, 4.923, 4.926, 4.929, 4.932, 4.935, 4.938, 4.941, 4.943, 4.946, 4.949, 4.951, 4.954, 4.956, 4.959, 4.961, 4.964, 4.966, 4.968, 4.971, 4.973, 4.975, 4.978, 4.980, 4.982, 4.984, 4.986, 4.988, 4.990, 4.992, 4.995, 4.996, 4.998, 5.000, 5.002, 5.004, 5.006, 5.008, 5.010, 5.011, 5.013, 5.014, 5.016, 5.017, 5.019, 5.020, 5.022, 5.023, 5.025, 5.026, 5.027, 5.028, 5.030, 5.031, 5.032, 5.033, 5.035, 5.036, 5.037, 5.038, 5.039, 5.040, 5.041, 5.043, 5.044, 5.045, 5.046, 5.047, 5.048, 5.049, 5.050, 5.051, 5.052, 5.053, 5.054, 5.055, 5.056, 5.057, 5.058, 5.059, 5.060, 5.061, 5.062, 5.062, 5.063, 5.064, 5.065, 5.066, 5.067, 5.068, 5.069, 5.070, 5.071, 5.072, 5.072, 5.073, 5.074, 5.075, 5.076, 5.077, 5.078, 5.079, 5.079, 5.080, 5.081, 5.082, 5.083, 5.084, 5.085, 5.085, 5.086, 5.087, 5.088, 5.089, 5.090, 5.090, 5.091, 5.092, 5.093, 5.094, 5.094, 5.095, 5.096, 5.097, 5.098, 5.098, 5.099, 5.100, 5.101, 5.101, 5.102, 5.103, 5.104, 5.104, 5.105, 5.106, 5.107, 5.107, 5.108, 5.109, 5.109, 5.110, 5.111, 5.111, 5.112, 5.113, 5.114, 5.114, 5.115, 5.116, 5.116, 5.117, 5.117, 5.118, 5.119, 5.119, 5.120, 5.121, 5.121, 5.122, 5.122, 5.123, 5.124, 5.124, 5.125, 5.125, 5.126, 5.127, 5.127, 5.128, 5.128, 5.129, 5.130, 5.130, 5.131, 5.131, 5.132, 5.132, 5.133, 5.133, 5.134, 5.135, 5.135, 5.136, 5.136, 5.137, 5.137, 5.138, 5.138, 5.139, 5.139, 5.140, 5.140, 5.141, 5.141, 5.142, 5.142, 5.143, 5.143, 5.144, 5.144, 5.145, 5.145, 5.146, 5.146, 5.147, 5.147, 5.148, 5.148, 5.149, 5.149, 5.150, 5.150, 5.151, 5.151, 5.152, 5.152, 5.152, 5.153, 5.153, 5.154, 5.154, 5.155, 5.155, 5.156, 5.156, 5.156, 5.157, 5.157, 5.158, 5.158, 5.159, 5.159, 5.160, 5.160, 5.160, 5.161, 5.161, 5.162, 5.162, 5.162, 5.163, 5.163, 5.164, 5.164, 5.164, 5.165, 5.165, 5.166, 5.166, 5.166, 5.167, 5.167, 5.168, 5.168, 5.168, 5.169, 5.169, 5.170, 5.170, 5.170, 5.171, 5.171, 5.171, 5.172, 5.172, 5.172, 5.173, 5.173, 5.174, 5.174, 5.174, 5.175, 5.175, 5.175, 5.176, 5.176, 5.176, 5.177, 5.177, 5.177, 5.178, 5.178, 5.179, 5.179, 5.179, 5.180, 5.180, 5.180, 5.181, 5.181, 5.181, 5.182, 5.182, 5.182, 5.183, 5.183, 5.183, 5.184, 5.184, 5.184, 5.184, 5.185, 5.185, 5.185, 5.186, 5.186, 5.186, 5.187, 5.187, 5.187, 5.188, 5.188, 5.188, 5.189, 5.189, 5.189, 5.189, 5.190, 5.190, 5.190, 5.191, 5.191, 5.191, 5.192, 5.192, 5.192, 5.192, 5.193, 5.193, 5.193, 5.194, 5.194, 5.194, 5.194, 5.195, 5.195, 5.195, 5.196, 5.196, 5.196, 5.196, 5.197, 5.197, 5.197, 5.198, 5.198, 5.198, 5.198, 5.199, 5.199, 5.199, 5.199, 5.200, 5.200, 5.200, 5.201, 5.201, 5.201, 5.201, 5.202, 5.202, 5.202, 5.202, 5.203, 5.203, 5.203, 5.203, 5.204, 5.204, 5.204, 5.204, 5.205, 5.205, 5.205, 5.205, 5.206, 5.206, 5.206, 5.206, 5.207, 5.207, 5.207, 5.207, 5.208, 5.208, 5.208, 5.208, 5.209, 5.209, 5.209, 5.209, 5.210, 5.210, 5.210, 5.210, 5.211, 5.211, 5.211, 5.211, 5.212, 5.212, 5.212, 5.212, 5.213, 5.213, 5.213, 5.213, 5.213, 5.214, 5.214, 5.214, 5.214, 5.215, 5.215, 5.215, 5.215, 5.215, 5.216, 5.216, 5.216, 5.216, 5.217, 5.217, 5.217, 5.217, 5.217, 5.218, 5.218, 5.218, 5.218, 5.219, 5.219, 5.219, 5.219, 5.219, 5.220, 5.220, 5.220, 5.220, 5.220, 5.221, 5.221, 5.221, 5.221, 5.222, 5.222, 5.222, 5.222, 5.222, 5.223, 5.223, 5.223, 5.223, 5.223, 5.224, 5.224, 5.224, 5.224, 5.224, 5.225, 5.225, 5.225, 5.225, 5.225, 5.226, 5.226, 5.226, 5.226, 5.226, 5.227, 5.227, 5.227, 5.227, 5.227, 5.228, 5.228, 5.228, 5.228, 5.228, 5.229, 5.229, 5.229, 5.229, 5.229, 5.230, 5.230, 5.230, 5.230, 5.230, 5.231, 5.231, 5.231, 5.231, 5.231, 5.232, 5.232, 5.232, 5.232, 5.232, 5.233, 5.233, 5.233, 5.233, 5.233, 5.234, 5.234, 5.234, 5.234, 5.234, 5.234, 5.235, 5.235, 5.235, 5.235, 5.235, 5.236, 5.236, 5.236, 5.236, 5.236, 5.237, 5.237, 5.237, 5.237, 5.237, 5.237, 5.238, 5.238, 5.238, 5.238, 5.238, 5.239, 5.239, 5.239, 5.239, 5.239, 5.239, 5.240, 5.240, 5.240, 5.240, 5.240, 5.241, 5.241, 5.241, 5.241, 5.241, 5.241, 5.242, 5.242, 5.242, 5.242, 5.242, 5.242, 5.243, 5.243, 5.243, 5.243, 5.243, 5.243, 5.244, 5.244, 5.244, 5.244, 5.244, 5.245, 5.245, 5.245, 5.245, 5.245, 5.245, 5.246, 5.246, 5.246, 5.246, 5.246, 5.246, 5.247, 5.247, 5.247, 5.247, 5.247, 5.247, 5.248, 5.248, 5.248, 5.248, 5.248, 5.248, 5.248, 5.249, 5.249, 5.249, 5.249, 5.249, 5.249, 5.250, 5.250, 5.250, 5.250, 5.250, 5.250, 5.251, 5.251, 5.251, 5.251, 5.251, 5.251, 5.252, 5.252, 5.252, 5.252, 5.252, 5.252, 5.252, 5.253, 5.253, 5.253, 5.253, 5.253, 5.253, 5.254, 5.254, 5.254, 5.254, 5.254, 5.254, 5.254, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255, 5.255 };

/* H ionizing luminosity corresponding to constant star formation at 1
   Msun/yr, expressed in photons / sec (derived from a slug
   calculation at an age of 50 Myr */
const double qH0perSFR = 1.57e53;

/* H ionizing luminosity, expressed as log_10 photons per second per 10^6 Msun worth of stars */
const double qH0_SB99[1000] =
  {  52.615, 52.615, 52.613, 52.613, 52.613, 52.613, 52.613, 52.614, 52.614, 52.615, 52.615, 52.616, 52.615, 52.614, 52.615, 52.614, 52.614, 52.614, 52.613, 52.612, 52.608, 52.600, 52.603, 52.605, 52.608, 52.609, 52.611, 52.612, 52.615, 52.612, 52.610, 52.596, 52.596, 52.600, 52.601, 52.600, 52.600, 52.600, 52.599, 52.591, 52.586, 52.579, 52.565, 52.550, 52.545, 52.537, 52.525, 52.509, 52.492, 52.467, 52.451, 52.433, 52.407, 52.390, 52.363, 52.341, 52.295, 52.270, 52.303, 52.342, 52.366, 52.385, 52.393, 52.403, 52.401, 52.385, 52.363, 52.351, 52.335, 52.317, 52.303, 52.289, 52.270, 52.247, 52.219, 52.199, 52.186, 52.141, 52.104, 52.089, 52.070, 52.047, 52.037, 52.019, 51.997, 51.984, 51.968, 51.943, 51.923, 51.903, 51.890, 51.865, 51.850, 51.815, 51.769, 51.750, 51.744, 51.720, 51.702, 51.676, 51.672, 51.657, 51.602, 51.570, 51.540, 51.524, 51.496, 51.470, 51.446, 51.422, 51.397, 51.376, 51.359, 51.349, 51.323, 51.285, 51.267, 51.245, 51.226, 51.203, 51.177, 51.150, 51.134, 51.118, 51.095, 51.078, 51.052, 51.031, 51.013, 50.990, 50.970, 50.951, 50.932, 50.911, 50.892, 50.874, 50.856, 50.837, 50.820, 50.800, 50.783, 50.767, 50.745, 50.727, 50.708, 50.690, 50.673, 50.655, 50.640, 50.622, 50.600, 50.583, 50.565, 50.556, 50.541, 50.520, 50.505, 50.491, 50.478, 50.457, 50.442, 50.428, 50.412, 50.395, 50.383, 50.366, 50.354, 50.340, 50.323, 50.312, 50.298, 50.282, 50.272, 50.257, 50.248, 50.237, 50.219, 50.208, 50.198, 50.186, 50.172, 50.160, 50.154, 50.140, 50.128, 50.112, 50.102, 50.087, 50.077, 50.063, 50.054, 50.038, 50.031, 50.019, 50.010, 49.999, 49.993, 49.978, 49.972, 49.956, 49.946, 49.938, 49.926, 49.918, 49.908, 49.899, 49.888, 49.875, 49.869, 49.855, 49.847, 49.835, 49.828, 49.816, 49.804, 49.797, 49.785, 49.770, 49.768, 49.755, 49.748, 49.736, 49.727, 49.721, 49.713, 49.701, 49.694, 49.684, 49.677, 49.666, 49.655, 49.649, 49.638, 49.630, 49.624, 49.613, 49.602, 49.596, 49.587, 49.576, 49.567, 49.560, 49.549, 49.541, 49.531, 49.524, 49.515, 49.510, 49.499, 49.490, 49.479, 49.477, 49.464, 49.453, 49.450, 49.443, 49.431, 49.427, 49.419, 49.411, 49.404, 49.402, 49.391, 49.383, 49.373, 49.370, 49.361, 49.353, 49.344, 49.341, 49.330, 49.322, 49.313, 49.310, 49.302, 49.290, 49.283, 49.280, 49.271, 49.259, 49.253, 49.245, 49.240, 49.233, 49.222, 49.220, 49.209, 49.200, 49.197, 49.188, 49.180, 49.177, 49.166, 49.158, 49.153, 49.149, 49.146, 49.138, 49.132, 49.125, 49.120, 49.114, 49.106, 49.100, 49.097, 49.090, 49.082, 49.077, 49.068, 49.063, 49.055, 49.052, 49.044, 49.039, 49.031, 49.025, 49.024, 49.012, 49.006, 49.002, 48.998, 48.988, 48.985, 48.978, 48.972, 48.967, 48.955, 48.950, 48.949, 48.941, 48.936, 48.929, 48.922, 48.916, 48.912, 48.908, 48.902, 48.896, 48.887, 48.882, 48.880, 48.874, 48.869, 48.860, 48.855, 48.849, 48.846, 48.846, 48.837, 48.837, 48.830, 48.823, 48.821, 48.814, 48.811, 48.804, 48.800, 48.795, 48.789, 48.787, 48.783, 48.778, 48.776, 48.762, 48.766, 48.758, 48.752, 48.749, 48.747, 48.736, 48.735, 48.731, 48.726, 48.720, 48.715, 48.714, 48.709, 48.705, 48.697, 48.695, 48.691, 48.682, 48.681, 48.675, 48.670, 48.664, 48.660, 48.655, 48.649, 48.645, 48.643, 48.638, 48.633, 48.629, 48.624, 48.618, 48.613, 48.607, 48.603, 48.603, 48.597, 48.591, 48.587, 48.583, 48.579, 48.576, 48.569, 48.567, 48.564, 48.557, 48.554, 48.549, 48.544, 48.541, 48.539, 48.533, 48.530, 48.526, 48.525, 48.523, 48.519, 48.513, 48.508, 48.504, 48.499, 48.495, 48.495, 48.491, 48.487, 48.482, 48.481, 48.473, 48.472, 48.464, 48.461, 48.465, 48.458, 48.454, 48.452, 48.445, 48.441, 48.437, 48.432, 48.427, 48.426, 48.419, 48.417, 48.419, 48.415, 48.409, 48.403, 48.399, 48.391, 48.389, 48.384, 48.386, 48.382, 48.377, 48.370, 48.367, 48.364, 48.361, 48.357, 48.353, 48.352, 48.342, 48.343, 48.337, 48.333, 48.329, 48.322, 48.323, 48.315, 48.313, 48.309, 48.305, 48.303, 48.298, 48.293, 48.289, 48.284, 48.284, 48.281, 48.280, 48.277, 48.273, 48.267, 48.270, 48.263, 48.260, 48.256, 48.250, 48.246, 48.244, 48.241, 48.239, 48.236, 48.238, 48.231, 48.229, 48.221, 48.219, 48.216, 48.212, 48.210, 48.206, 48.201, 48.197, 48.198, 48.194, 48.188, 48.190, 48.184, 48.178, 48.175, 48.174, 48.170, 48.164, 48.161, 48.155, 48.157, 48.156, 48.148, 48.151, 48.141, 48.136, 48.133, 48.129, 48.126, 48.122, 48.120, 48.113, 48.113, 48.109, 48.106, 48.103, 48.098, 48.095, 48.091, 48.087, 48.084, 48.076, 48.073, 48.071, 48.068, 48.065, 48.060, 48.056, 48.059, 48.050, 48.047, 48.042, 48.039, 48.033, 48.030, 48.024, 48.024, 48.020, 48.017, 48.011, 48.008, 48.005, 48.003, 48.001, 47.999, 47.993, 47.991, 47.988, 47.989, 47.985, 47.981, 47.978, 47.976, 47.972, 47.968, 47.966, 47.960, 47.954, 47.952, 47.949, 47.944, 47.948, 47.941, 47.938, 47.934, 47.932, 47.928, 47.924, 47.917, 47.914, 47.912, 47.912, 47.906, 47.903, 47.902, 47.898, 47.894, 47.891, 47.888, 47.886, 47.884, 47.878, 47.876, 47.869, 47.867, 47.864, 47.861, 47.857, 47.854, 47.853, 47.847, 47.844, 47.845, 47.841, 47.842, 47.831, 47.827, 47.822, 47.819, 47.812, 47.812, 47.809, 47.806, 47.801, 47.799, 47.797, 47.794, 47.788, 47.786, 47.783, 47.778, 47.774, 47.771, 47.768, 47.761, 47.762, 47.759, 47.753, 47.751, 47.751, 47.746, 47.744, 47.740, 47.734, 47.731, 47.729, 47.723, 47.722, 47.715, 47.712, 47.710, 47.704, 47.707, 47.701, 47.703, 47.695, 47.689, 47.685, 47.685, 47.679, 47.678, 47.680, 47.670, 47.671, 47.668, 47.665, 47.663, 47.662, 47.660, 47.656, 47.651, 47.652, 47.650, 47.647, 47.643, 47.637, 47.636, 47.632, 47.634, 47.626, 47.624, 47.619, 47.618, 47.611, 47.614, 47.613, 47.614, 47.600, 47.602, 47.597, 47.594, 47.592, 47.589, 47.586, 47.588, 47.576, 47.577, 47.572, 47.572, 47.568, 47.566, 47.571, 47.557, 47.556, 47.552, 47.551, 47.546, 47.551, 47.544, 47.540, 47.533, 47.527, 47.528, 47.525, 47.523, 47.516, 47.514, 47.508, 47.508, 47.508, 47.503, 47.500, 47.495, 47.499, 47.495, 47.490, 47.494, 47.486, 47.477, 47.474, 47.470, 47.468, 47.465, 47.462, 47.457, 47.452, 47.449, 47.450, 47.444, 47.444, 47.442, 47.437, 47.437, 47.434, 47.428, 47.426, 47.424, 47.415, 47.416, 47.411, 47.408, 47.402, 47.403, 47.399, 47.389, 47.391, 47.388, 47.377, 47.380, 47.380, 47.377, 47.373, 47.370, 47.367, 47.364, 47.361, 47.359, 47.356, 47.351, 47.344, 47.344, 47.340, 47.338, 47.332, 47.333, 47.331, 47.320, 47.320, 47.319, 47.317, 47.308, 47.306, 47.306, 47.302, 47.295, 47.298, 47.294, 47.292, 47.294, 47.292, 47.284, 47.286, 47.283, 47.280, 47.273, 47.275, 47.280, 47.271, 47.272, 47.264, 47.263, 47.265, 47.263, 47.264, 47.258, 47.253, 47.248, 47.250, 47.242, 47.239, 47.237, 47.230, 47.232, 47.229, 47.222, 47.223, 47.218, 47.225, 47.216, 47.213, 47.210, 47.208, 47.205, 47.199, 47.201, 47.195, 47.193, 47.193, 47.198, 47.188, 47.181, 47.181, 47.178, 47.174, 47.171, 47.171, 47.169, 47.166, 47.168, 47.171, 47.169, 47.154, 47.149, 47.146, 47.143, 47.140, 47.133, 47.134, 47.131, 47.136, 47.122, 47.127, 47.125, 47.120, 47.113, 47.111, 47.111, 47.106, 47.103, 47.107, 47.104, 47.100, 47.100, 47.093, 47.090, 47.087, 47.086, 47.090, 47.081, 47.074, 47.072, 47.075, 47.072, 47.078, 47.065, 47.057, 47.067, 47.065, 47.056, 47.049, 47.044, 47.041, 47.039, 47.039, 47.036, 47.034, 47.028, 47.026, 47.025, 47.022, 47.016, 47.016, 47.016, 47.010, 47.009, 47.006, 47.004, 47.002, 46.997, 46.987, 46.982, 46.979, 46.977, 46.977, 46.974, 46.972, 46.974, 46.972, 46.965, 46.963, 46.964, 46.961, 46.960, 46.952, 46.952, 46.945, 46.943, 46.939, 46.938, 46.935, 46.933, 46.934, 46.929, 46.923, 46.922, 46.923, 46.918, 46.916, 46.913, 46.911, 46.911, 46.905, 46.903, 46.898, 46.894, 46.891, 46.888, 46.887, 46.885, 46.883, 46.876, 46.872, 46.872, 46.866, 46.864, 46.862, 46.860, 46.858, 46.854, 46.850, 46.845, 46.843, 46.840, 46.836, 46.829, 46.827, 46.826, 46.829, 46.825, 46.826, 46.824, 46.824, 46.822, 46.819, 46.831, 46.828, 46.814, 46.812, 46.808, 46.807, 46.807, 46.799, 46.800, 46.800, 46.799, 46.800, 46.795, 46.793, 46.787, 46.786, 46.788, 46.786, 46.792, 46.781, 46.777, 46.775, 46.773, 46.770, 46.772, 46.764, 46.758, 46.760, 46.758, 46.756, 46.749, 46.756, 46.750, 46.772, 46.760, 46.745 };

void
userEOS(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	void *params,
	double *gamma, double *delta) {
  fprintf(stderr,
	  "Warning: userEOS function called but not implemented!\n");
  return;
}

void
userAlpha(const double t, const double dt, const grid *grd,
	  const double *col, const double *pres, const double *eInt,
	  const double *gamma, const double *delta,
	  void *params,
	  double *alpha) {

  gsl_poly_complex_workspace *wksp;
  double coef[6], sol[10];
  double omega, kappa, T1, sigma, kcrit, J, Q, d, tfast, tgrowth;
  unsigned long i, j;
  const int m=2;
  gsl_complex nu2, nu;
  double alphacoef = ((paramType *) params)->alphacoef;

  /* Allocate workspace */
  wksp = gsl_poly_complex_workspace_alloc(6);

  /* Store the coefficiencts of the dispersion relation polynomial
     that don't change */
  coef[2] = -8.0;
  coef[3] = coef[4] = 0.0;

  /* Loop over cells */
  for (i=0; i<grd->nr; i++) {

    /* Handle negative column densities, which may arise if the time
       step overshoots. We just need to make sure the code doesn't
       barf here, as the negative column densities will be fixed by
       the iterative solver. */
    if (col[i] <= 0.0) {
      alpha[i] = ALPHAMIN;
      continue;
    }

    /* The various quantities needed for the stability analysis */
    sigma = sqrt(pres[i]/col[i]);
    omega = grd->vphi_g[i+1]/grd->r_g[i+1];
    kappa = sqrt(2.0*(1.0+grd->beta_g[i+1]))*omega;
    T1 = -(2.0*m*omega/(kappa*grd->r_g[i+1])) *
      (2.0*m*omega/(kappa*grd->r_g[i+1])) * (grd->beta_g[i+1]-1);
    kcrit = kappa*kappa / (2.0*M_PI*G*col[i]);
    Q = kappa*sigma / (M_PI*G*col[i]);
    if (Q < 0)
      printf("Q < 0! kappa = %f, sigma = %f, col = %f\n", kappa, sigma, col[i]);
    J = sqrt(T1)/kcrit;

    /* Coefficients of the disperison relation polynomial */
    coef[0] = -Q*Q*Q*Q;
    coef[1] = 6.0*Q*Q;
    coef[5] = 16.0*J*J;

    /* Find roots, giving minima and maxima of D; note that we test
       to make sure that the leading coefficient isn't too small,
       which makes the problem ill-conditioned as causes the solver
       to hang. If we detect this happening, we just drop the leading
       coefficient and solve a 4th order polynomial instead of a 5th. */
    if (coef[5] > 1e-30) {
      gsl_poly_complex_solve(coef, 6, wksp, sol);
    } else {
      gsl_poly_complex_solve(coef, 5, wksp, sol);
    }

    /* For each root, compute the growth time of the instability */
    tfast = 1.0e30;
    for (j=0; j<5; j++) {

      /* Skip roots with negative real part, or imaginary part that is
	 greater than roundoff */
      if ((sol[2*j] < 0) || (fabs(sol[2*j+1]) > 1.0e-6)) continue;

      /* Compute D */
      d = (Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j]) *
	(Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j] -
	 4.0*J*J*sol[2*j]*sol[2*j]);

      /* Ignore stable roots */
      if (d > 0) continue;

      /* Get frequency and growth time */
      GSL_SET_COMPLEX(&nu2,
		      1.0 + 0.5*(Q*Q/(4.0*sol[2*j]*sol[2*j]) - 1.0/sol[2*j]),
		      0.5 * sqrt(-d));
      nu = gsl_complex_sqrt(nu2);
      tgrowth = 1.0 / (GSL_IMAG(nu)*kappa) / (2.0*M_PI/omega);

      /* Store minimum growth time */
      if (tgrowth < tfast) tfast = tgrowth;
    }

    /* Now check for gravitational instability */
    if (Q < 1) {
      GSL_SET_COMPLEX(&nu, 0.0, sqrt(1.0/(Q*Q)-1.0));
      tgrowth = 1.0 / (GSL_IMAG(nu)*kappa) / (2.0*M_PI/omega);
      if (tgrowth < tfast) tfast = tgrowth;
    }

    /* Compute alpha based on the fastest growing mode timescale */
    alpha[i] = alphacoef * exp(-tfast+1.0);
    if (alpha[i] > 1.0) alpha[i] = 1.0;
    if (alpha[i] < ALPHAMIN) alpha[i] = ALPHAMIN;
  }

  /* Free workspace */
  gsl_poly_complex_workspace_free(wksp);
}

void
userMassSrc(const double t, const double dt, const grid *grd,
	    const double *col, const double *pres, const double *eInt,
	    const double *gamma, const double *delta,
	    void *params,
	    double *massSrc) {
  /* Estimate Mdot, including effects of star formation and mass ejection */

  long i;
  double *sigma_sfr, *sigma_wind;
  double col_floor = ((paramType *) params)->col_floor;

  /* Allocate memory to hold results */
  sigma_sfr = calloc(grd->nr, sizeof(double));
  sigma_wind = calloc(grd->nr, sizeof(double));

  /* Get SFR and wind mass flux */
  sfrwind(grd, col, pres, params, sigma_sfr, sigma_wind);

  /* Compute mass source term */
  for (i=0; i<grd->nr; i++) massSrc[i] = -sigma_sfr[i] - sigma_wind[i];

  /* Apply absolute floor */
  if (col_floor > 0) {
    for (i=0; i<grd->nr; i++) {
      massSrc[i] += col_floor / (grd->r_g[i+1]/grd->vphi_g[i+1]) *
	exp(col_floor/col[i]) / (1.0 + exp(SQR(col[i]/col_floor)));
    }
  }

  /* Free memory */
  free(sigma_sfr);
  free(sigma_wind);
}

void
userIntEnSrc(const double t, const double dt, const grid *grd,
	     const double *col, const double *pres, const double *eInt,
	     const double *gamma, const double *delta,
	     void *params, double *intEnSrc) {
  /* Cooling rate = eta Sigma sigma^2 Omega = eta P vphi/r; heating
     rate from star formation = pdot * sigma / area */
  long i;
  double eta = ((paramType *) params)->eta;
  double sigmath = ((paramType *) params)->sigmath;
  double shapefac = ((paramType *) params)->shapefac;
  double zetad = ((paramType *) params)->zetad;
  double sigma_max = ((paramType *) params)->sigma_max;
  double col_min = ((paramType *) params)->col_min;
  double *pdot = PDOT(params, grd);
  double sigma2, sigmaNT, rhostar, a, b, c, h;

  for (i=0; i<grd->nr; i++) {
    /* Gas velocity dispersion */
    sigma2 = pres[i]/col[i];
    if (sigma2 > SQR(sigmath)) {
      sigmaNT = sqrt(sigma2 - SQR(sigmath));
      /* Stellar density */
      rhostar = shapefac * SQR(grd->vphi_g[i+1]) *
	(1.0+2.0*grd->beta_g[i+1]) /
	(4.0*M_PI*G*SQR(grd->r_g[i+1]));
      /* Get scale height from OML model */
      if (rhostar > 0.0) {
	a = 2.0*M_PI*zetad*G*rhostar*col[i];
	b = M_PI/2.0*G*SQR(col[i]);
	c = -pres[i];
	h = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
      } else {
	h = pres[i] / (M_PI/2.0*G*SQR(col[i]));
      }
      /* Cooling rate = eta Sigma sigmaNT^2 / (h/sigmaNT) */
      intEnSrc[i] = -eta * col[i] * SQR(sigmaNT) / (h/sigmaNT);
    } else {
      intEnSrc[i] = 0.0;
    }

    /* Heating rate = pdot/area * sigma */
    intEnSrc[i] += pdot[i] * sqrt(sigma2) * exp(-col_min/col[i]);

    /* Extra cooling in hot cells */
    intEnSrc[i] -= col[i] * SQR(sigmaNT) /
      (grd->r_g[i+1]/grd->vphi_g[i+1]) *
      exp(sigmaNT/sigma_max) / (1.0 + exp(SQR(sigma_max/sigmaNT)));
  }
}

void
userIBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type ibc_pres, const enth_bc_type ibc_enth,
	void *params,
	double *ibc_pres_val, double *ibc_enth_val) {
  fprintf(stderr,
	  "Warning: userIBC function called but not implemented!\n");
  return;
}

void
userOBC(const double t, const double dt, const grid *grd,
	const double *col, const double *pres, const double *eInt,
	const double *gamma, const double *delta,
	const pres_bc_type obc_pres, const enth_bc_type obc_enth,
	void *params,
	double *obc_pres_val, double *obc_enth_val) {
  fprintf(stderr,
	  "Warning: userOBC function called but not implemented!\n");
  return;
}

void
userPreTimestep(const double t, const double dt,
		const grid *grd, double *col, double *pres,
		double *eInt, double *mBnd, double *eBnd,
		double *mSrc, double *eSrc,
		void *params, const unsigned long nUserOut,
		double *userOut) {
  /* We do two things before a time step: (1) record the stellar mass
     at the start of the time step so we can back out the mass lost
     during it; (2) compute and store the current momentum and energy
     injection rate from the stellar population */
  long i, j;
  double p_per_sn = ((paramType *) params)->p_per_sn;
  double e_per_sn = ((paramType *) params)->e_per_sn;
  double star_rad_var = ((paramType *) params)->star_rad_var;
  double *msrc_tmp = MSRC_TMP(params);
  double *pdot = PDOT(params, grd);
  double *edot = EDOT(params, grd);
  double *mdot = MDOT(params, grd);
  double *pdot_tmp, *pdot_tmp1, *edot_tmp, *edot_tmp1,
    *mdot_tmp, *mdot_tmp1;
  double pdot_tot, norm;
  double *Sigmastar = SIGMASTAR(params, grd);

  /* Record cumulative mass source term up to this point */
  for (i=0; i<grd->nr; i++) msrc_tmp[i] = mSrc[i];

  /* Compute pdot and edot from star formation */
  pdot_tmp = calloc(grd->nr, sizeof(double));
  edot_tmp = calloc(grd->nr, sizeof(double));
  mdot_tmp = calloc(grd->nr, sizeof(double));
  pdot_tmp1 = calloc(grd->nr, sizeof(double));
  edot_tmp1 = calloc(grd->nr, sizeof(double));
  mdot_tmp1 = calloc(grd->nr, sizeof(double));
  for (i=0; i<grd->nr; i++)
    injection_rate(Sigmastar+i*nbin_SB99, p_per_sn, e_per_sn,
		   pdot_tmp+i, edot_tmp+i, mdot_tmp+i);

  /* Initialize to 0 */
  for (i=0; i<grd->nr; i++) pdot[i] = edot[i] = mdot[i] = 0.0;

  /* Loop over cells to do convolution, being careful to ensure that
     the area-integrated momentum and energy added stays constant */
  for (j=0; j<grd->nr; j++) {
    pdot_tot = 0.0;
    for (i=0; i<grd->nr; i++) {
      pdot_tmp1[i] = exp(-SQR(grd->r_g[i+1]-grd->r_g[j+1]) /
			 (2.0*SQR(star_rad_var*grd->r_g[j+1])))
	* pdot_tmp[j];
      pdot_tot += pdot_tmp1[i] * grd->area[i];
      edot_tmp1[i] = exp(-SQR(grd->r_g[i+1]-grd->r_g[j+1]) /
			 (2.0*SQR(star_rad_var*grd->r_g[j+1])))
	* edot_tmp[j];
      mdot_tmp1[i] = exp(-SQR(grd->r_g[i+1]-grd->r_g[j+1]) /
			 (2.0*SQR(star_rad_var*grd->r_g[j+1])))
	* mdot_tmp[j];
    }
    norm = pdot_tmp[j]*grd->area[j] / (pdot_tot + SMALL);
    for (i=0; i<grd->nr; i++) {
      pdot[i] += norm * pdot_tmp1[i];
      edot[i] += norm * edot_tmp1[i];
      mdot[i] += norm * mdot_tmp1[i];
    }
  }

  /* Free memory */
  free(pdot_tmp);
  free(pdot_tmp1);
  free(edot_tmp);
  free(edot_tmp1);
  free(mdot_tmp);
  free(mdot_tmp1);
}

void
userPostTimestep(const double t, const double dt,
		 const grid *grd, double *col, double *pres,
		 double *eInt, double *mBnd, double *eBnd,
		 double *mSrc, double *eSrc,
		 void *params, const unsigned long nUserOut,
		 double *userOut) {
  long i, j;
  long *tctr = &(((paramType *) params)->tctr);
  double *msrc_tmp = MSRC_TMP(params);
  double *Sigmastar = SIGMASTAR(params, grd);
  double *pdot = PDOT(params, grd);
  double *edot = EDOT(params, grd);
  double *mdot = MDOT(params, grd);
  double *Sigmastar_tot = userOut;
  double *Sigmawind_tot = userOut + grd->nr;
  double *pdot_save = userOut + 2*grd->nr;
  double *edot_save = userOut + 3*grd->nr;
  double *sigma_sfr = userOut + 4*grd->nr;
  double *sigma_wind = userOut + 5*grd->nr;
  double *sigma_sfr_obs = userOut + 6*grd->nr;
  double *mdot_save = userOut + 7*grd->nr;
  double d2qdtdA;

  /* Safety check */
  assert(nUserOut == 8);

  /* Shift bins of star formation history if needed */
  while (t/(dt_SB99*MYR) >= *tctr+1) {
    for (i=0; i<grd->nr; i++) {
      for (j=nbin_SB99-1; j>0; j--) {
	Sigmastar[j+i*nbin_SB99] = Sigmastar[j-1+i*nbin_SB99];
      }
      Sigmastar[i*nbin_SB99] = 0.0;
    }
    (*tctr)++;
  }

  /* Figure out the current SFR and mass loss rate; we will use these
     to partition the change in mass between new stars and wind */
  sfrwind(grd, col, pres, params, sigma_sfr, sigma_wind);

  /* Record stars formed during this timestep in the youngest bin;
     also add to source terms we're tracking */
  for (i=0; i<grd->nr; i++) {
    Sigmastar[i*nbin_SB99] += (-mSrc[i] + msrc_tmp[i]) *
      sigma_sfr[i] / (sigma_sfr[i] + sigma_wind[i] + SMALL);
  }

  /* Save cumulative surface densities of stars formed and wind mass
     lost, and instantaneous momentum injection rate */
  for (i=0; i<grd->nr; i++) {
    Sigmastar_tot[i] += (-mSrc[i] + msrc_tmp[i]) *
      sigma_sfr[i] / (sigma_sfr[i] + sigma_wind[i] + SMALL);
    Sigmawind_tot[i] += (-mSrc[i] + msrc_tmp[i]) *
      sigma_wind[i] / (sigma_sfr[i] + sigma_wind[i] + SMALL);
    pdot_save[i] = pdot[i];
    edot_save[i] = edot[i];
    mdot_save[i] = mdot[i];
  }

  /* Save "observed" SFR at this time, i.e., SFR that would be
     inferred from Halpha or similar ionization-based tracer */
  for (i=0; i<grd->nr; i++) {
    d2qdtdA = 0.0;
    for (j=0; j<nbin_SB99; j++)
      d2qdtdA += pow(10.0, qH0_SB99[j]) *
	(Sigmastar[i*nbin_SB99+j] / (1.0e6*MSUN));
    sigma_sfr_obs[i] = d2qdtdA / (qH0perSFR/(MSUN/YR));
  }
}

/* This is an auxiliary routine to compute the momentum and energy
   input per unit area from a stellar population */
void injection_rate(double *Sigmastar, double p_per_sn,
		    double e_per_sn, double *pdot_tot,
		    double *edot_tot, double *mdot_tot) {
  long i;
  double snr = 0.0;
  double pdot_sn = 0.0, pdot_wind = 0.0, pdot_rad = 0.0;
  double edot_sn = 0.0, edot_wind = 0.0;

  /* Figure out the total SN rate and momentum injection */
  for (i=0; i<nbin_SB99; i++) {
    if ((snr_SB99[i] > -29.0) && (Sigmastar[i] > 0.0)) {
      snr += pow(10.0, snr_SB99[i]) * Sigmastar[i]/(1.0e6*MSUN)/YR;
    }
  }
  pdot_sn = snr*p_per_sn;
  edot_sn = snr*e_per_sn;

  /* Get wind momentum and energy input rate */
  for (i=0; i<nbin_SB99; i++) {
    pdot_wind += pow(10.0, pdot_wind_SB99[i]) *
      (Sigmastar[i] / (1.0e6*MSUN));
    edot_wind += pow(10.0, edot_wind_SB99[i]) *
      (Sigmastar[i] / (1.0e6*MSUN));
  }

  /* Get radiation momentum flux */
  for (i=0; i<nbin_SB99; i++)
    pdot_rad += pow(10.0, lum_SB99[i]) * (Sigmastar[i] / (1.0e6*MSUN)) / C;

  /* Set momentum and energy total values */
  *pdot_tot = pdot_rad + pdot_sn + pdot_wind;
  *edot_tot = edot_sn + edot_wind;

  /* Get mass injection rate */
  *mdot_tot = pow(10.0, mass_return_SB99[0]) / (dt_SB99*MYR) *
      (Sigmastar[0] / 1.0e6);
  for (i=1; i<nbin_SB99; i++)
    *mdot_tot += (pow(10.0, mass_return_SB99[i]) -
		  pow(10.0, mass_return_SB99[i-1])) / (dt_SB99*MYR) *
      (Sigmastar[i] / 1.0e6);
}

/* This routine takes the current state of the simulation as input and
   computes the star formation rate and wind mass loss rate in every cell */
void sfrwind(const grid *grd, const double *col, const double *pres,
	     void *params, double *sigma_SFR, double *sigma_wind) {
  long i;
  double shapefac = ((paramType *) params)->shapefac;
  double zetad = ((paramType *) params)->zetad;
  double epsff0 = ((paramType *) params)->epsff0;
  double sigmath = ((paramType *) params)->sigmath;
  double col_min = ((paramType *) params)->col_min;
  double en_wind = ((paramType *) params)->en_wind;
  double *pdot = PDOT(params, grd);
  double *edot = EDOT(params, grd);
  double *mdot = MDOT(params, grd);
  double rhostar, a, b, c, h;
  double rhog, tff, alphavir, epsff;
  double pdotEdd, xcrit, sigma2, mach2, r, sigma_surf2, zeta, vwind;

  /* Loop over cells */
  for (i=0; i<grd->nr; i++) {

    /* Stellar density */
    rhostar = shapefac * SQR(grd->vphi_g[i+1]) *
      (1.0+2.0*grd->beta_g[i+1]) /
      (4.0*M_PI*G*SQR(grd->r_g[i+1]));

    /* Get scale height from OML model */
    if (rhostar > 0.0) {
      a = 2.0*M_PI*zetad*G*rhostar*col[i];
      b = M_PI/2.0*G*SQR(col[i]);
      c = -pres[i];
      h = (-b + sqrt(b*b-4.0*a*c))/(2.0*a);
    } else {
      h = pres[i] / (M_PI/2.0*G*SQR(col[i]));
    }

    /* Gas density and free-fall time */
    rhog = col[i]/(2*h);
    tff = sqrt(3*M_PI/(32*G*rhog));

    /* Virial ratio */
    alphavir = pres[i]/(M_PI/2.0*G*col[i]*col[i]*h);
    /*alphavir = tff/(h/sqrt(pres[i]/col[i]));*/

    /* epsff */
    epsff = exp(alphavir*log(epsff0));
    /*epsff = exp(sqrt(alphavir)*log(epsff0));*/

    /* SFR */
    sigma_SFR[i] = epsff * col[i] / tff;

    /* Compute wind using either momentum or energy prescription */
    if (en_wind <= 0.0) {

      /* Momentum case */

      /* Eddington momentum injection rate */
      /*pdotEdd = 2.0*M_PI*G*col[i]*(col[i]+rhostar*h);*/
      pdotEdd = 2.0*M_PI*G*col[i]*(col[i]+rhostar*grd->r_g[i+1]/shapefac);

      /* xcrit */
      xcrit = log(pdot[i]/pdotEdd);

      /* Dispersion of surface densities */
      sigma2 = pres[i]/col[i];
      mach2 = sigma2 / SQR(sigmath);
      r = 0.5 * (3.0-2.5)/(2.0-2.5) * (1.0-pow(mach2, 2.0-2.5)) /
	(1.0-pow(mach2, 3.0-2.5));
      sigma_surf2 = log(1.0 + r*mach2/4.0);

      /* zeta */
      zeta = 0.5*(1.0-erf((-2.0*xcrit+sigma_surf2) / sqrt(8.0*sigma_surf2)));

      /* Wind mass flux per unit area,  which we take to be zeta *
	 Sigma_gas / t_dyn, where t_dyn = h / sigma */
      sigma_wind[i] = zeta * col[i] * sqrt(sigma2)/h * exp(-col_min/col[i]);

    } else {

      /* Energy case */

      /* Thermal velocity of hot gas, assuming complete thermalization */
      vwind = sqrt(2.0*edot[i]/(mdot[i]+SMALL));

      /* Wind mass launch rate, using Martin 2005 entrainment model */
      sigma_wind[i] = sqrt(3.0*mdot[i]*rhog*vwind);

      /*if ((i==285) && (sigma_wind[i] > 1e-16))
	printf("sigma_sfr = %e, sigma_wind = %e\n", sigma_wind[i], sigma_SFR[i]);*/

    }
  }
}

void
userCheckRead(
	      FILE *fp, grid *grd, const unsigned long nOut,
	      double *tOut, double *colOut,
	      double *presOut, double *eIntOut, double *mBndOut,
	      double *eBndOut, double *mSrcOut, double *eSrcOut,
	      const unsigned long nUserOut, double *userOut,
	      void *params
	      ) {
  long *tctr = &(((paramType *) params)->tctr);
  double *Sigmastar = SIGMASTAR(params, grd);

  /* Read counter and the age-dependent stellar surface densities */
  fread(tctr, sizeof(long), 1, fp);
  fread(Sigmastar, nbin_SB99*grd->nr, sizeof(double), fp);
}

void
userCheckWrite(
	      FILE *fp,
	      const grid *grd, const unsigned long nOut,
	      const double *tOut, const double *colOut,
	      const double *presOut, const double *eIntOut,
	      const double *mBndOut, const double *eBndOut,
	      const double *mSrcOut, const double *eSrcOut,
	      const unsigned long nUserOut, const double *userOut,
	      const void *params
	      ) {
  long *tctr = &(((paramType *) params)->tctr);
  double *Sigmastar = SIGMASTAR(params, grd);

  /* Write the counter and the age-dependent stellar surface densities */
  fwrite(tctr, sizeof(long), 1, fp);
  fwrite(Sigmastar, nbin_SB99*grd->nr, sizeof(double), fp);
}

max_stars_repo_name: franciscaconcha/amuse-vader

max_issues_repo_name: franciscaconcha/amuse-vader

max_forks_repo_name: franciscaconcha/amuse-vader


text:
#ifndef QDM_THETA_H
#define QDM_THETA_H 1

#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>

int
qdm_theta_optimize(
    gsl_vector *result,
    gsl_vector *emperical_quantiles,
    gsl_matrix *ix
);

void
qdm_theta_matrix_constrain(
    gsl_matrix *theta,
    double min
);

#endif /* QDM_THETA_H */

max_stars_repo_name: calebcase/qdm

max_issues_repo_name: calebcase/qdm

max_forks_repo_name: calebcase/qdm


text:
/***********************************************************************
* mexeigPartial.c : C mex file for the routine DSYEVX in LAPACK
*
* mex -O -largeArrayDims -lmwlapack -lmwblas mexeigPartial.c
*
*  [V,D,nv] = mexeigPartial(A, vl, vu, options);
*
* options.jobz = 1 (default), compute both eigenvectors and eigenvalues;
*              = 2, compute only eigenvalues.
* options.range = 1 (default), compute all;
*               = 2, compute  eigenvalues between [vl, vu].
*               = 3, compute  vl-th to vu-th eigenvalues.
* options.abstol = tolerance (default = 1e-8)
***********************************************************************/

#include <math.h>
#include <mex.h>
#include <matrix.h>
#include <blas.h>
#include <lapack.h>
#include <string.h> /* needed for memcpy() */

#if !defined(_WIN32)
#define dsyevd dsyevd_
#define dsyevx dsyevx_
#define dgesdd dgesdd_
#endif

/**********************************************************
*
***********************************************************/
void mexFunction(
      int nlhs,   mxArray  *plhs[],
      int nrhs,   const mxArray  *prhs[] )

{    double   *A, *AA, *nV, *V, *D, VL, VU, abstol, *work, *work2, *tmp;

     ptrdiff_t   m, n, IL, IU, lwork, *iwork, *ifail, info, options;
     char     *jobz="V";
     char     *uplo="U";
     char     *range="A";

/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

   if (nrhs > 4){
      mexErrMsgTxt("mexeigPartial: requires at most 4 input arguments."); }
   if (nlhs > 3){
      mexErrMsgTxt("mexeigPartial: requires at most 3 output argument."); }

/* CHECK THE DIMENSIONS */

    m = mxGetM(prhs[0]);
    n = mxGetN(prhs[0]);
    if (m != n) {
       mexErrMsgTxt("mexeigPartial: matrix must be square."); }
    if (mxIsSparse(prhs[0])) {
       mexErrMsgTxt("mexeigPartial: sparse matrix not allowed."); }
    A = mxGetPr(prhs[0]);

    options = 1;
    abstol = 1e-8;
    /*Get fields from options */
    tmp = mxGetPr(mxGetField(prhs[3], 0, "abstol"));  abstol  = *tmp;
    tmp = mxGetPr(mxGetField(prhs[3], 0, "jobz"));    options = (ptrdiff_t)*tmp;
    tmp = mxGetPr(mxGetField(prhs[3], 0, "range"));   options = (ptrdiff_t)*tmp;
    if (options==1) {
       range="A";  }
    else if (options==2) {
       range="V"; VL = *mxGetPr(prhs[1]); VU = *mxGetPr(prhs[2]);  }
    else if (options == 3) {
       range="I";
       IL = (ptrdiff_t)*mxGetPr(prhs[1]);
       IU = (ptrdiff_t)*mxGetPr(prhs[2]);
    }
    /***** create return argument *****/
    plhs[0] = mxCreateDoubleMatrix(n,n,mxREAL);
    V = mxGetPr(plhs[0]);
    plhs[1] = mxCreateDoubleMatrix(n,1,mxREAL);
    D = mxGetPr(plhs[1]);
    plhs[2] = mxCreateDoubleMatrix(1,1,mxREAL);
    nV = mxGetPr(plhs[2]);

    /***** Do the computations in a subroutine *****/
    lwork  = 8*n;
    work   =  mxCalloc(lwork,sizeof(double));
    iwork  =  mxCalloc(5*n, sizeof(ptrdiff_t));
    ifail  =  mxCalloc(n, sizeof(ptrdiff_t));

    AA =  mxCalloc(n*n,sizeof(double));
    memcpy(AA,mxGetPr(prhs[0]),(n*n)*sizeof(double));

    dsyevx(jobz, range, uplo, &n, AA, &n, &VL, &VU, &IL, &IU,
           &abstol, &m, D, V, &n, work, &lwork, iwork, ifail, &info);

    *nV = m;
    return;
 }
/**********************************************************/

max_stars_repo_name: napinoco/BBCPOP

max_issues_repo_name: napinoco/BBCPOP

max_forks_repo_name: napinoco/BBCPOP


text:
/**
 */
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <gsl/gsl_interp.h>
#include <gsl/gsl_multifit_nlin.h>
#include "aXe_grism.h"
#include "aXe_utils.h"
//#include "spce_PET.h"
#include "inout_aper.h"
//#include "trace_conf.h"
#include "aper_conf.h"
//#include "spc_sex.h"
//#include "disp_conf.h"
//#include "spc_wl_calib.h"
//#include "spce_binning.h"
#include "spc_spc.h"
//#include "spce_output.h"
//#include "spc_FITScards.h"
#include "fringe_conf.h"
#include "trfit_utils.h"
#include "ipixcorr_utils.h"

#define AXE_IMAGE_PATH "AXE_IMAGE_PATH"
#define AXE_OUTPUT_PATH "AXE_OUTPUT_PATH"
#define AXE_CONFIG_PATH "AXE_CONFIG_PATH"

int
main (int argc, char *argv[])
{
  char *opt;

  char ipc_aper_file[MAXCHAR];
  char ipc_aper_file_path[MAXCHAR];

  char aper_file[MAXCHAR];
  char aper_file_path[MAXCHAR];

  char conf_file[MAXCHAR];
  char conf_file_path[MAXCHAR];

  char grism_file[MAXCHAR];
  char grism_file_path[MAXCHAR];

  char bck_image[MAXCHAR];
  char bck_image_path[MAXCHAR];

  char IPC_file[MAXCHAR];
  char IPC_file_path[MAXCHAR];

  char ipc_file[MAXCHAR];
  char ipc_file_path[MAXCHAR];

  char SPC_file[MAXCHAR];
  char SPC_file_path[MAXCHAR];

  char IPC_func[MAXCHAR];
  char IPC_func_path[MAXCHAR];

  aperture_conf *conf;

  object **oblist;
  observation *obs=NULL, *bck;

  //Initialize obs struct memebers
  obs->grism = NULL;
  obs->pixerrs = NULL;
  obs->dq = NULL;

  interpolator *ipcorr;

  interpolator *nlincorr;

  int index;
  //int i;

  full_spectr *SPC;

  fitsfile  *SPC_ptr;
  fitsfile  *IPC_ptr;
  FITScards *cards;

  int f_status=0;

  int aperID=0, beamID=0;
  //objindex;

  int spec_OAF=0;
  int point_like=0;

  beam act_beam;
  beam *beam_ptr;

  double max_ext=0.0;
  //double adcgain;
  double exptime=NAN;

  gsl_matrix *data_matrix;

  int ipixc=0;
  int nlinc=0;
  int icorr=0;
  int num=0;

  if ((argc < 3) || (opt = get_online_option ("help", argc, argv)))
    {
      fprintf (stdout,
         "aXe_FRINGECORR Version %s:\n",RELEASE);
      exit (1);
    }

  fprintf (stdout, "aXe_INTPIXCORR: Starting...\n");

  // reset the parameter index
  index = 0;

  // read in the grism image name
  strcpy(grism_file, argv[++index]);
  build_path (AXE_IMAGE_PATH, grism_file, grism_file_path);

  // read in the configuration file name
  strcpy(conf_file, argv[++index]);
  build_path (AXE_CONFIG_PATH, conf_file, conf_file_path);

  // load the configuration file
  conf = get_aperture_descriptor(conf_file_path);

  // determine the science extension number
  get_extension_numbers(grism_file_path, conf,conf->optkey1,conf->optval1);

  // get or set up the input OAF name
  if ((opt = get_online_option("in_OAF", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(aper_file, opt);

      // mark that there is a special OAF file
      spec_OAF=1;
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(aper_file, grism_file);
      replace_file_extension (grism_file, aper_file, ".fits",
			      ".OAF", conf->science_numext);
    }
  build_path(AXE_OUTPUT_PATH, aper_file, aper_file_path);

  // build the name of the modified
  // OAF file
  replace_file_extension (grism_file, ipc_aper_file, ".fits",
			      "_ipc.OAF", conf->science_numext);
  build_path(AXE_OUTPUT_PATH, ipc_aper_file, ipc_aper_file_path);

  // get or set up the output SPC name
  if ((opt = get_online_option("out_SPC", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(IPC_file, opt);
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(IPC_file, grism_file);
      replace_file_extension(grism_file, IPC_file, ".fits",
			     "_ipc.SPC.fits", conf->science_numext);
    }
  build_path(AXE_OUTPUT_PATH, IPC_file, IPC_file_path);

  // get or set up the output SPC name
  if ((opt = get_online_option("in_SPC", argc, argv)))
    {
      // copy the parameter value if given
      strcpy(SPC_file, opt);
    }
  else
    {
      // copy the grism name and construct the name
      // if not given as parameter
      strcpy(SPC_file, grism_file);
      replace_file_extension(grism_file, SPC_file, ".fits",
			     ".SPC.fits", conf->science_numext);
    }
  build_path(AXE_OUTPUT_PATH, SPC_file, SPC_file_path);

  // form the background image name
  replace_file_extension(grism_file, bck_image, ".fits",
			 ".BCK.fits", conf->science_numext);
  build_path (AXE_OUTPUT_PATH, bck_image, bck_image_path);

  // check whether an intra pixel sensitivity correction
  // shall be applied
  if ((opt = get_online_option("ipixcorr", argc, argv)))
    ipixc = 1;
  // check whether an intra pixel sensitivity correction
  // shall be applied
  if ((opt = get_online_option("nlincorr", argc, argv)))
    nlinc = 1;

  // check whether something is
  // done at all
  if (ipixc || nlinc)
    {
      // Build the intrapixel correction file name
      sprintf(IPC_func,"%s",conf->IPIXfunc);

      // overwrite it with the parameter if given
      if ((opt = get_online_option ("IPIXFUNCTION", argc, argv)))
	strcpy(IPC_func,opt);

      if (!strcmp(IPC_file,"None"))
	// check whether the correction function is given;
	// complain if not given!
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_INTPIXCORR: Correction function must be given!!\n");
      else
	// compose the full pathname
	build_path (AXE_CONFIG_PATH, IPC_func, IPC_func_path);

      // get the maximum extension number
      if ((opt = get_online_option ("max_ext", argc, argv)))
	max_ext = atof(opt);

      // check whether either max_ext is defined
      // OR a special OAF file is given
      if (!max_ext && !spec_OAF)
	aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		     "aXe_INTPIXCORR: Either maximum extension\n"
		     "OR a special OAF file must be given!\n");
    }

  // print on screen what's done
  fprintf (stdout, "aXe_INTPIXCORR: Grism file name:                     %s\n",
	   grism_file_path);
  fprintf (stdout, "aXe_INTPIXCORR: Configuration file name:             %s\n",
	   conf_file_path);
  fprintf (stdout, "aXe_INTPIXCORR: Correction function:                 %s\n",
	   IPC_func_path);
  fprintf (stdout, "aXe_INTPIXCORR: Input OAF file name:                 %s\n",
	   aper_file_path);
  fprintf (stdout, "aXe_INTPIXCORR: Input SPC file name:                 %s\n",
	   SPC_file_path);
  fprintf (stdout, "aXe_INTPIXCORR: Output SPC file name:                %s\n",
	   IPC_file_path);
  fprintf (stdout, "aXe_INTPIXCORR: Maximal extension for point objects: %f\n",
	   max_ext);
  if (ipixc)
    fprintf (stdout, "aXe_INTPIXCORR: Intra-pixel sensitivity correction is aplied!\n");
  if (nlinc)
    fprintf (stdout, "aXe_INTPIXCORR: Non-linearity correction is aplied!\n");

  // check whether something is
  // done at all
  if (ipixc || nlinc)
    {
      // check whether the intra pixel
      // sensitivity is applied
      if (ipixc)
	{
	  // determine the exposure time
	  if (isnan(exptime))
	    exptime = (double)get_float_from_keyword(grism_file_path, 1,
						     conf->exptimekey);
	  if (isnan(exptime))
	    exptime = 1.0;

	  // load the image data
	  obs = load_image_t(grism_file_path, conf->science_numext,
			     conf->errors_numext, conf->dq_numext,
			     conf->dqmask, exptime, conf->rdnoise);

	  // load the background data
	  bck = load_image_t(bck_image_path, conf->science_numext,
			     conf->errors_numext, conf->dq_numext,
			     conf->dqmask, exptime, conf->rdnoise);

	  // subtract the background from
	  // the image data
	  data_matrix = bcksub_observation(obs, bck);
	}

      // load the aperture list
      fprintf (stdout, "aXe_INTPIXCORR: Loading object aperture list...");fflush(stdout);
      oblist = file_to_object_list_seq (aper_file_path, obs);
      fprintf (stdout,"%d objects loaded.\n",object_list_size(oblist));fflush(stdout);

      // open the input SPC file
      SPC_ptr = get_SPC_opened(SPC_file_path, 0);

      // open the output SPC file
      IPC_ptr = create_SPC_opened (IPC_file_path,1);
      cards = get_FITS_cards_opened(SPC_ptr);
      put_FITS_cards_opened(IPC_ptr, cards);
      free_FITScards(cards);

      if (ipixc)
	// create an interpolator reading in the fits file
	ipcorr = create_interp_ftable(IPC_func_path,2,"X","FX",IPCORR_INTERP_TYPE);
      else
	ipcorr = NULL;

      if (nlinc)
	// for the NICMOS nonlinearity correction
	// create an interpolator for the non-linear correction
	//nlincorr = create_nlincor();
	nlincorr = create_interp_ftable(IPC_func_path,2,"LAMBDA","BPAR",NLINCORR_INTERP_TYPE);
      else
	nlincorr = NULL;

      while ((aperID!=-1) || (beamID!=-1))
	{
	  // load the next extension from the SPC table
	  SPC = get_ALL_from_next_in_SPC(SPC_ptr, &aperID, &beamID);

	  // break the loop if the end was reached
	  if ((aperID==-1) && (beamID==-1))
	    break;

	  // if the full spectrum is NULL,
	  // go to the next beam
	  if (!SPC)
	    continue;

	  // select the actual beam from the beam list
	  act_beam = find_beam_in_object_list(oblist, aperID, beamID);
	  beam_ptr = find_beamptr_in_object_list(oblist, aperID, beamID);

	  // determine whether the object is pointlike
	  point_like = is_pointlike(act_beam, spec_OAF, max_ext);

	  // if the obejct is indeed
	  // a pointlike one
	  if (point_like)
	    {
	      // check whether the intra-pixel sensititvity correction
	      // shall be applied
	      if (ipixc)
		{
		  // form the name of the file which contain
		  // the data
		  replace_file_extension(grism_file, ipc_file, ".fits",
					 ".IPC.dat", aperID);
		  build_path(AXE_OUTPUT_PATH, ipc_file, ipc_file_path);

		  icorr = fitting_ipc_corr(act_beam, conf_file_path, ipcorr,
					   SPC, obs, data_matrix, ipc_file_path, beam_ptr);

		  if (icorr)
		    fprintf(stdout, "aXe_INTPIXCORR: Corrected BEAM %d%c\n", SPC->aperID, BEAM(SPC->beamID));
		  else
		    fprintf(stdout, "aXe_INTPIXCORR: phase accuracy in BEAM %d%c too low\n", SPC->aperID, BEAM(SPC->beamID));

		}

	      // check whether the non-linearity correction
	      // shall be applied
	      if (nlinc)
		{
		  /*
		  if (!strcmp(conf->gainkey,"None"))
		    adcgain=1.0;
		  else
		    adcgain = (double)get_float_from_keyword(grism_file_path, 1, conf->gainkey);
		  */
		  // apply the NICMOS non-linearity correction
		  fprintf(stdout, "aXe_INTPIXCORR: Non-linearity correting BEAM %d%c...\n", SPC->aperID, BEAM(SPC->beamID));
		  // the gain mostly used in NICMOS
		  // is hardcoded here!!!!!
		  nlin_corr_beam(nlincorr, 6.5, SPC);
		}
	    }
	  else
	    {
	      fprintf(stdout, "aXe_INTPIXCORR: Skipping BEAM %d%c...\n", SPC->aperID, BEAM(SPC->beamID));
	    }

	  add_spectra_to_SPC_opened (IPC_ptr, SPC->fgr_spec, SPC->bck_spec,
				     SPC->obj_spec, SPC->aperID, SPC->beamID);

	  /* Copy header from OPET extension into this SPC extension */
	  cards = get_FITS_cards_opened(SPC_ptr);
	  put_FITS_cards_opened(IPC_ptr,cards);
	  free_FITScards(cards);

	  free_full_spectr(SPC);
      	}

      refurbish_object_list(oblist, 1, -100, 0);

      fprintf (stdout, "aXe_INTPIXCORR: Writing aperture file %s ... ", ipc_aper_file);fflush(stdout);
      num = object_list_to_file (oblist, ipc_aper_file_path, 0);
      fprintf (stdout, "%d beams written.\n", num);

      // free the object list
      if (oblist!=NULL)
	free_oblist (oblist);

      if (ipcorr)
	// free the interpolator
	free_interp(ipcorr);

      if (nlincorr)
	// free the interpolator
	free_interp(nlincorr);

      fits_close_file (SPC_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_INTPIXCORR: " "Error closing SPC: %s\n",
		       SPC_file_path);
	}

      fits_close_file (IPC_ptr, &f_status);
      if (f_status)
	{
	  ffrprt (stderr, f_status);
	  aXe_message (aXe_M_FATAL, __FILE__, __LINE__,
		       "aXe_INTPIXCORR: " "Error closing SPC: %s\n",
		       IPC_file_path);
	}
    }
  // say goodbye
  fprintf (stdout, "aXe_INTPIXCORR: Done...\n");
  exit (0);
}

max_stars_repo_name: sosey/pyaxe

max_issues_repo_name: sosey/pyaxe

max_forks_repo_name: sosey/pyaxe