#!/usr/bin/perl
#
# Make a simple model of an absorption line in a stellar spectrum.
#   Model the line as a gaussian.
#
# Print to stdout lines that look like this:
#
#         wavelength      intensity
#
#
# MWR 9/10/2008
#
# This is modified to compute some quantities which might tell us
#    whether gaussian or lorentzian is a better fit.
#
# MWR 9/18/2008

$debug = 0;

$C = 2.9979E8;
if ($debug > 0) {
  printf " speed of light C = %13.8e \n", $C;
}

# read in all the data
$infile = "./spec_closeup.dat";
open(INFILE, "$infile") || die("can't open file $infile");
$num_data = 0;
while (<INFILE>) {
  $line = " " . $_;
  chomp($line);
  @words = split(/\s+/, $line);
  $data_lam[$num_data] = $words[1];
  $data_inten[$num_data] = $words[2];
  if ($debug > 1) {
    printf "  read line %6d  lam %9.3f  inten %8.1f \n", 
               $num_data, $data_lam[$num_data], $data_inten[$num_data];
  }

  $num_data++;
}
close(INFILE);


$start_lambda = 6545.0;
$end_lambda = 6547.0;
$d_lambda = 0.01;

# these are the parameters of the model
#    intensity units are arbitrary linear 
$continuum_level = 9680.0;
#    wavelength units are Angstroms
$min_line_intensity = 6730.0;
$central_lambda = 6546.2475;
$fwhm_lambda = 0.135;
$sigma = $fwhm_lambda / 2.35;

# compute the central frequency (Hz)
$central_freq = ($central_lambda*1.0E-10);
if ($debug > 0) {
  printf " central_freq is %13.8e \n", $central_freq;
}
$central_freq = $C / ($central_lambda*1.0E-10);
if ($debug > 0) {
  printf " central_freq is %13.8e \n", $central_freq;
}

# assume that the uncertainty in each datum is 0.5 units of intensity
$uncert_inten = 0.5;

# for each of the wavelengths with a measurement, compute
#      the value of our model.

$start_chisq_lam = 6546.12;
$end_chisq_lam = 6546.38;
$num_chisq_points = 0;

$chisq = 0.0;
for ($i = 0; $i < $num_data; $i++) {

  if (($data_lam[$i] < $start_chisq_lam) || ($data_lam[$i] > $end_chisq_lam)) {
    next;
  }

  # get the frequency of this datum
  $freq = $C / ($data_lam[$i] * 1.0E-10);
  if ($debug > 0) {
    printf "     this_freq is %13.8e \n", $freq;
  }

  if ($debug > 0) {
    printf "  lam %13.8e  cont %13.8e  inten %13.8e \n",
             $data_lam[$i], $continuum_level, $data_inten[$i];
  }

  # here we test for the Lorentzian
  $lor_left_side = ($continuum_level) / ($continuum_level - $data_inten[$i]);
  $lor_right_side = ($freq - $central_freq)*($freq - $central_freq);

  # now here we test for gaussian
  $gauss_left_side = ($continuum_level - $data_inten[$i]) / ($continuum_level);
  $gauss_left_side = log($gauss_left_side);
  $gauss_right_side = ($freq - $central_freq)*($freq - $central_freq);


  printf "%13.8e %13.8e %13.8e %13.8e %13.8e   \n", 
        $data_lam[$i], $data_inten[$i], $lor_right_side, $lor_left_side,
        $gauss_left_side;

}
#printf "#  chisq is %9.4e based on %6d data\n", $chisq, $num_chisq_points;



exit 0;



##########################################################################
# PROCEDURE: compute_model
#
# DESCRIPTION: Given a wavelength and four parameters, compute the intensity 
#            of a model absorption line at the given wavelength.
#
#            Call with arguments like so:
#
#                compute_model(lam, C, D, lam0, fwhm)
#            
#            The wavelength 'lam' should be in Angstroms.
#
#            The four parameters are:
#
#                 C         level of continuum (arb linear intensity)
#                 D         max depth of line (arb linear intensity)
#                 lam0      central wavelength of line (Angstroms)
#                 fwhm      FWHM of line (Angstroms)
#
#            The model is
#
#                                       -  (lam - lam0)^2 / 2*sigma^2
#                 model  =  C  -  D * e
#
#            where
#                       sigma  =  FWHM / 2.35
#
# RETURNS:
#            value of the model at the given wavelength
#
sub compute_model {

  my($C, $D, $lam0, $fwhm, $lam);
  my($x, $arg, $normalized_depth, $model_intensity);
  my($sigma);

  $lam = $_[0];
  $C = $_[1];
  $D = $_[2];
  $lam0 = $_[3];
  $fwhm = $_[4];

  $sigma = $fwhm / 2.35;

  $x = $lam - $lam0;
  $arg = ($x*$x) / (2.0*$sigma*$sigma);
  $normalized_depth = exp(0.0 - $arg);
  $model_intensity = $C - $D*$normalized_depth;

  return($model_intensity);
}



