/*//////////////////////////////////////////////////////////////////////////////
// <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 species mass fraction at a given time
//       in multiple GCE models (zones).
//     </abstract>
//   </description>
// </file>
//////////////////////////////////////////////////////////////////////////////*/
 
#include <WnSimpleGce.h>

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

  WnSimpleGce **p_zones;
  WnSimpleGce__Species *p_species;
  double d_initial_gas_mass, d_Delta, d_tmp = 0., d_omega, d_primary_yield;
  unsigned int i, i_k;
  size_t i_zones, i_yield_coefficients;
  FILE *p_file;
  char s_species[32];

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

   if( argc != 3 ) {
      fprintf(
        stderr,
        "\nUsage: %s file_name \n\n",
        argv[0]
      );
      fprintf(
        stderr,
        "file_name = name of the input file\n\n"
      );
      fprintf(
        stderr,
        "t = time in Gyr at which to compute mass fractions\n\n"
      );
      return EXIT_FAILURE;
   }

  /*============================================================================
  // Open file.
  //==========================================================================*/

  p_file = fopen( argv[1], "r" );

  if( !p_file )
  {
    fprintf( stderr, "Couldn't open file.\n" );
    return EXIT_FAILURE;
  }

  /*============================================================================
  // Create species.
  //==========================================================================*/

  fscanf( p_file, "%s\n", s_species );

  p_species = WnSimpleGce__Species__new( s_species );

  fscanf( p_file, "%lu\n", (unsigned long *) &i_yield_coefficients );

  for( i = 0; i < i_yield_coefficients; i++ )
  {
    fscanf( p_file, "%lf\n", &d_tmp );
    WnSimpleGce__Species__updateYieldCoefficient( p_species, i, d_tmp );
  }

  /*============================================================================
  // Read and update the decay rate.
  //==========================================================================*/

  fscanf( p_file, "%lf\n", &d_tmp );

  WnSimpleGce__Species__updateDecayRate( p_species, d_tmp );

  /*============================================================================
  // Read the primary yield.
  //==========================================================================*/

  fscanf( p_file, "%lf\n", &d_primary_yield );

  /*============================================================================
  // Create zones.
  //==========================================================================*/

  fscanf( p_file, "%lu\n", (unsigned long *) &i_zones );

  p_zones = 
    ( WnSimpleGce ** ) malloc( sizeof( WnSimpleGce * ) * i_zones );

  /*============================================================================
  // Read in zone data.
  //==========================================================================*/

  for( i = 0; i < i_zones; i++ )
  {
    fscanf(
      p_file,
      "%lf %lf %lf %ud\n",
      &d_initial_gas_mass,
      &d_Delta,
      &d_omega,
      &i_k
    );

    p_zones[i] = WnSimpleGce__new( );

    WnSimpleGce__updateInitialGasMass( p_zones[i], d_initial_gas_mass );
  
    WnSimpleGce__updateInfallKValue( p_zones[i], i_k );
  
    WnSimpleGce__updateInfallDelta( p_zones[i], d_Delta );
  
    WnSimpleGce__updateOmega( p_zones[i], d_omega );

    WnSimpleGce__updatePrimaryMetallicityYield( p_zones[i], d_primary_yield );

  }

  fclose( p_file );
  
  /*============================================================================
  // Compute species mass fraction.
  //==========================================================================*/

  fprintf(
    stdout,
    "\nt (Gyr) = %12.6e\n",
    atof( argv[2] )
  );

  for( i = 0; i < i_zones; i++ )
    fprintf(
      stdout,
      "Zone = %d  X(%s) = %.4e\n",
      i,
      WnSimpleGce__Species__getName( p_species ),
      WnSimpleGce__computeSpeciesMassFraction(
        p_zones[i],
        p_species,
        atof( argv[2] )
      )
    );

  fprintf( stdout, "\n" );

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

  WnSimpleGce__Species__free( p_species );

  for( i = 0; i < i_zones; i++ )
    WnSimpleGce__free( p_zones[i] );

  free( p_zones );

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

  return EXIT_SUCCESS;

}