/*//////////////////////////////////////////////////////////////////////////////
// <file type="public">
//   <license>
//     This file was originally written by Bradley S. Meyer.
//
//     This is free software; you can redistribute it and/or modify it
//     under the terms of the GNU General Public License as published by
//     the Free Software Foundation; either version 2 of the License, or
//     (at your option) any later version.
//
//     This software 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.
//
//     Please see the src/README.txt file in this distribution for more
//     information.
//   </license>
//   <description>
//     <abstract>
//       Example to compute the remainder of a radioactive species
//       (the mass fraction of the species in the gas relative to that
//       if the species were stable) in a simple GCE model from input
//       yield and lifetime data.
//     </abstract>
//   </description>
// </file>
//////////////////////////////////////////////////////////////////////////////*/
 
#include <WnSimpleGce.h>

#define D_INTERVAL 0.5

int
main( int argc, char **argv ) {

  WnSimpleGce *p_model;
  WnSimpleGce__Species *p_rad, *p_rad_as_stable;
  double d_t, d_rad, d_stable;

  /*============================================================================
  // Check input.
  //==========================================================================*/

   if ( argc != 9 ) {
      fprintf(
        stderr,
        "\nUsage: %s k Delta omega alpha species decay_rate alpha beta\n\n",
        argv[0]
      );
      fprintf(
        stderr,
        "k = Clayton infall parameter k\n\n"
      );
      fprintf(
        stderr,
        "Delta = Clayton infall parameter Delta\n\n"
      );
      fprintf(
        stderr,
        "omega = gas consumption rate (per Gyr)\n\n"
      );
      fprintf(
        stderr,
        "alpha = primary metallicity yield\n\n"
      );
      fprintf(
        stderr,
        "species = name of radioactive species\n\n"
      );
      fprintf(
        stderr,
        "decay_rate = decay rate (per Gyr) of radioactive species\n\n"
      );
      fprintf(
        stderr,
        "alpha = alpha parameter for species\n\n"
      );
      fprintf(
        stderr,
        "beta = beta parameter for species\n\n"
      );
      return EXIT_FAILURE;
   }

  /*============================================================================
  // Create model.
  //==========================================================================*/

  p_model = WnSimpleGce__new( );

  /*============================================================================
  // Update model.
  //==========================================================================*/

  WnSimpleGce__updateInfallKValue( p_model, (unsigned int) atoi( argv[1] ) );
  
  WnSimpleGce__updateInfallDelta( p_model, atof( argv[2] ) );
  
  WnSimpleGce__updateOmega( p_model, atof( argv[3] ) );
  
  WnSimpleGce__updatePrimaryMetallicityYield( p_model, atof( argv[4] ) );

  /*============================================================================
  // Create species and stable version and assign and yield coefficients.
  //==========================================================================*/

  p_rad = WnSimpleGce__Species__new( argv[5] );

  WnSimpleGce__Species__updateYieldCoefficient( p_rad, 0, atof( argv[7] ) );

  WnSimpleGce__Species__updateYieldCoefficient( p_rad, 1, atof( argv[8] ) );

  p_rad_as_stable = WnSimpleGce__Species__new( argv[5] );

  WnSimpleGce__Species__updateYieldCoefficient(
    p_rad_as_stable,
    0,
    atof( argv[7] )
  );

  WnSimpleGce__Species__updateYieldCoefficient(
    p_rad_as_stable,
    1,
    atof( argv[8] )
  );

  /*============================================================================
  // Assign decay rate.
  //==========================================================================*/

  WnSimpleGce__Species__updateDecayRate( p_rad, atof( argv[6] ) );

  /*============================================================================
  // Compute mass fraction ratio.
  //==========================================================================*/

  fprintf(
    stdout,
    "  t (Gyr)\tr(%s)\n\n",
    WnSimpleGce__Species__getName( p_rad )
  );

  d_t = D_INTERVAL;

  while( d_t < 20.0 + D_INTERVAL ) {

    d_rad = WnSimpleGce__computeSpeciesMassFraction( p_model, p_rad, d_t );

    d_stable =
      WnSimpleGce__computeSpeciesMassFraction(
        p_model,
        p_rad_as_stable,
        d_t
    );

    fprintf(
      stdout,
      "%12.6e %12.6e\n",
      d_t,
      d_rad / d_stable
    );

    d_t += D_INTERVAL;

  }

  /*============================================================================
  // Do cleanup.
  //==========================================================================*/

  WnSimpleGce__Species__free( p_rad );
  WnSimpleGce__Species__free( p_rad_as_stable );
  WnSimpleGce__free( p_model );

  /*============================================================================
  // Done.
  //==========================================================================*/

  return EXIT_SUCCESS;

}