#include <stdio.h>
#include <stdarg.h>
#include <stdlib.h>
#include <string.h>
#include <strings.h>
#include <sys/types.h>
#include <fcntl.h>
#include <math.h>

/* 
Simulation of the motion of an expanding object within uniform ISM
==================================================================

This program is usable if the expansion of the object and its motion can
be treated independently (which implies that internal (explosion) energy
is large relative to the kinetic energy of the motion). Energy radiated
away (e.g. due to shock ionization) is not taken into account.

Math
----

The mass sweep-up rate (caused by the motion, not the expansion) can be
described as

    m_m'(t) = \pi \rho r(t) v(t) 

where $r(t)$ is the radius at time $t$ and $\rho$ is the mass density.
The velocity $v(t)$ of the object is obtained by the assumption, that
the the kinetic energy is preserved:

    v(t)^2 ( m_m(t) + m_0 ) = v(0)^2 m_0

where m_0 is the initial mass of the object. These two equations can be
combined into a single differential equation, which is solved
numerically by this program.

Usage
-----

Parameters are set by macros, see comments below. 

The expansion function `er(t)` can be re-defined and currently
implements the SNR evolution according to Blondin et al. 
(1998ApJ...500..342B). For the radiative phase, the declaration
parameter $\frac{1}{3}$  was used. The evolution model is an analytic
approximation of a 1D-simulation.

To use it, compile the program and run it:
    
    gcc ghostsim.c -lm -o ghostsim && ./ghostsim 
    
These columns are output
    t   : time in <time units>
    v   : motion velocity (not the expansion velocity) in 
          <length unit> / <time unit>
    s   : traveled distance in <length unit>
    m   : mass swept-up by motion (not by expansion) in <mass unit>
    r	: expansion radius
    m_e : mass swept-up by expansion in <mass unit>
    v_e : expansion velocity in <length unit> / <time unit>
*/

// implementation of the expansion function
// physical parameters
#define NA 1.7E-2	// ambient ISM density in atoms per cm^3
#define E51 1		// explosion energy normalized to 10^51 ergs

// transition time according to Blondin et al. in yr (1998ApJ...500..342B)
static double t_tr = 2.9e4 * exp( log(E51)*4.0/17.0 ) * exp( log(NA)*(-9.0/17.0) );
// Sedov radius at transition time in pc / yr
static double r_tr = -1;
// Sedov velocity at transition time in pc / yr
static double v_tr = -1;

// expansion radius in <length units> at time t in <time units>
static double er(double t) {	// t in years
    static double a = -1, b = -1;
    if ( t <= t_tr ) return 0.314 * exp( log(E51/NA*t*t) * 0.2 );	
    if ( b < 0 ) {
	b = 0.1256 * exp( log(E51/(NA*t_tr*t_tr*t_tr) ) * 0.2 ); // Sedov velocity at transition time in pc / yr
	printf("v_tr=%-10.4g\n",b);
	a = er(t_tr) - b*t_tr*3;
	b *= exp( log(t_tr) * 2/3 ) * 3;
    }
    return a + b*exp( log(t) / 3 );
}


// simulation parameters
// output step size in <time unit>
#define DT_O 3420

// number of output steps
#define O_N 100

// simulation steps per output step
#define DT_N 10000


// physical parameters
// density in <mass unit> per <length unit>^3
#define RHO (NA /*cm^-3*/ * 2.939e55 /*cm^3 pc^-3*/ * 1.661e-24 /*g*/ / 1.988e33 /*<g / <sun mass>*/ ) 	// in solarmass pc^-3

// define initial speed in <length unit> per <time unit>
#define V0 2.159e-4	// in pc / year

// define initial mass in <mass units>
#define M0 13.4	// in solarmass



void main() {
    double mm = 0;	// mass swept-up mass in <mass unit>
    double t = 0;	// time in <time units>
    double dt = (double)DT_O / (double)DT_N;
    double v = V0;	// velocity in <length unit> / <time units>
    double r;		// radius in <length unit>
    double s = 0;	// moved distance
    for (int i=0; i<O_N; i++) {
	for (int j=0; j<DT_N; j++) {
	    // pseudo implicit step using two fixpoint iterations for better stability
	    t += dt;
	    r = er(t);
	    double mm1 = mm;
	    for ( int k=0; k<2; k++ ) {
		v = (double)V0 * sqrt( (double)M0 / ( (double)M0 + mm1) );
		mm1 = mm + M_PI*r*r*v*(double)RHO*dt;
	    }
	    mm = mm1;
	    s += v*dt;
	}
	printf("t=%-10.4g v=%-10.4g s=%-10.4g m=%-10.4g  r=%-10.4g m_e=%-10.4g v_e=%-10.4g\n", t, v, s, mm,  r, 4.0/3.0*M_PI*r*r*r*RHO, (r-er(0.999*t))/(0.001*t));
    }
    printf("v0=%-10.4g m0=%-10.4g rho=%-10.4g\n", (double)V0, (double)M0, (double)RHO);
}