forked from tardis-sn/tardis
-
Notifications
You must be signed in to change notification settings - Fork 0
sample code
wkerzendorf edited this page Dec 8, 2011
·
2 revisions
##sl.c
/* This is a simple code to model the formation of a single spectral line
in a wind */
#include "sl.h"
int
main()
{
FILE *out_file;
int i, npackets2, n, n_spec;
double mid_freq, runtot, norm, norm2;
double bb(), zrand;
double this_e, this_freq, this_r, this_mu, this_nucmf;
double d_inner, d_outer, d_line;
double r_sob, tau_sob, tau_rand, x, nh_sob;
int this_dead;
double check1;
int nn;
double phi;
double dist_0;
npackets = 1000000;
lambda_min = 1000;
lambda_max = 1400;
phi = 1.e-13;
r0 = 6.96e10;
r1 = 10*6.96e10;
v1 = 0.01 * CLIGHT;
freq_0 = CLIGHT / 1216.10 / 1e-8;
freq_max = CLIGHT / lambda_min / 1e-8;
freq_min = CLIGHT / lambda_max / 1e-8;
dist_0 = 3.086e21;
/* start by setting up the randon number generator */
zseed = 28469387;
rng = gsl_rng_alloc (gsl_rng_ran3);
gsl_rng_set ( rng, zseed);
/* call it a few times to get it in motion. */
for (n = 0; n< 100; n++)
{
x = gsl_rng_uniform(rng);
}
runtot = 0;
for (i = 0; i < MBINS; i++)
{
spectrum.lower_freq[i] = freq_min + (i * (freq_max-freq_min) / (MBINS));
spectrum.delta_freq[i] =(freq_max-freq_min) / (MBINS);
mid_freq = spectrum.lower_freq[i] + (0.5*spectrum.delta_freq[i]);
spectrum.flux[i] = bb(mid_freq);
runtot = runtot + (spectrum.flux[i]*spectrum.delta_freq[i]);
}
printf("runtot %g, lum %g\n", runtot, 4.*PI*r0 * r0 * 5.6705e-5 * pow(40000., 4));
norm = npackets / runtot;
for (i = 0; i < MBINS; i++);
npackets2=0;
for (i = 0; i < MBINS; i++)
{
spectrum.nstart[i] = spectrum.flux[i]*spectrum.delta_freq[i] * norm;
npackets2 += spectrum.nstart[i];
spectrum.flux[i] = 0.0;
if (spectrum.nstart[i] < 1)
{
printf("Negative or 0 packets in input spectrum. %d. Abort.", spectrum.nstart[i]);
exit(0);
}
}
printf("Number of packets intended: %d\n", npackets);
printf("Number of packets used: %d\n", npackets2);
npackets = npackets2;
i = 0;
nn=0;
for (n = 0; n < npackets; n++)
{
if (nn < spectrum.nstart[i])
{
zrand = gsl_rng_uniform(rng);
this_freq = spectrum.lower_freq[i] + (zrand * spectrum.delta_freq[i]);
nn++;
}
else
{
nn=0;
i++;
zrand = gsl_rng_uniform(rng);
this_freq = spectrum.lower_freq[i] + (zrand * spectrum.delta_freq[i]);
nn++;
}
this_r = r0;
zrand = gsl_rng_uniform(rng);
this_mu = sqrt(zrand);
this_e = 1. / npackets2;
this_dead = 0;
while (this_dead < 1)
{
/* Compute distance from current point to edge of wind */
d_outer = (-1.*this_r * this_mu) + sqrt(r1*r1 + (this_r*this_r*((this_mu*this_mu) -1.)));
/* Compute distance from current point to hit core (if possible) */
check1 = r0*r0 + (this_r*this_r*((this_mu*this_mu) -1.));
if (check1 < 0)
{
/* no intersection with core */
d_inner = 1.e50;
}
else
{
if (this_mu < 0)
{
/* intersection is possible */
d_inner = (-1.*this_r * this_mu) - sqrt(check1);
}
else
{
/* intersection is in negative direction */
d_inner = 1.e50;
}
}
/* Compute distance to line */
this_nucmf = this_freq * (1. - (this_mu * v1 * this_r / r1 / CLIGHT));
if (this_nucmf < freq_0)
{
d_line = 1.1e50;
}
else
{
d_line = (this_nucmf - freq_0) / this_freq * CLIGHT * r1 / v1;
}
/* all the d's are now known so decide which is the relevant one */
if (d_outer < d_line)
{
/* It will leave without hitting the line */
if (d_outer < d_inner)
{
this_dead = 1;
/* This flags an escaping packet */
}
else
{
this_dead = 2;
/* this flags a packet which hit the core */
}
}
else
{
/* It has a chance to hit the line before escaping. */
if (d_inner < d_line)
{
/*It hit the core */
this_dead = 2;
}
else
{
/*It has a chance to hit the line */
/* Compute the t_sob of the line */
r_sob = sqrt( (this_r*this_r) + (d_line*d_line) + (2*this_r*d_line*this_mu));
nh_sob = phi /4 /PI/r_sob/r_sob/r_sob/v1*r1/MH*MSUN/YEAR;
tau_sob = 0.02655103 * 0.416 * nh_sob * CLIGHT * r1 / this_freq / v1;
/*Get a random tau */
zrand = gsl_rng_uniform(rng);
tau_rand = -1.*log(zrand);
if (tau_rand > tau_sob)
{
/* No line event */
if (d_inner > d_outer)
{
this_dead = 1;
}
else
{
this_dead = 2;
}
}
else
{
/*event happens */
/* Move and scatter packet */
this_r = r_sob;
this_e = this_e * (1. - (this_mu * this_r * v1 / r1 / CLIGHT));
zrand = gsl_rng_uniform(rng);
this_mu = (2.*zrand) - 1.;
this_freq = freq_0 / (1. - (this_mu * this_r * v1 / r1 / CLIGHT));
this_e = this_e / (1. - (this_mu * this_r * v1 / r1 / CLIGHT));
//this_dead = 2;
}
}
}
}
if (this_dead < 1 || this_dead >2)
{
printf("This_dead error. Abort\n");
exit(0);
}
if (this_dead == 1)
{
/*These are the escaping packets */
n_spec = (this_freq - spectrum.lower_freq[0]) / spectrum.delta_freq[0];
spectrum.flux[n_spec] += this_e;
}
}
if ((out_file = fopen("spec.out", "w+")) == NULL){
printf("Cannot open out_file.txt.\n");
exit(0);
}
for (i = 0; i < MBINS; i++)
{
norm2 = runtot /4. / PI / dist_0 / dist_0 / spectrum.delta_freq[i];
fprintf(out_file, "%g %g\n", spectrum.lower_freq[i] + 0.5*spectrum.delta_freq[i], spectrum.flux[i]*norm2);
}
fclose(out_file);
return 0;
}
/* The black-body function */
double
bb(nu)
double nu;
{
double temp;
double hnukt;
double ans;
temp = 40000;
hnukt = H * nu / KB / temp;
ans = 2. * H * nu * nu * nu / (exp(hnukt) - 1.) / CLIGHT / CLIGHT;
ans = ans * 4 * PI * r0 * r0 * PI;
return(ans);
}
##sl.h
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#define MPACKETS 1000000
#define MBINS 1000
#define CLIGHT 2.99792458e10
#define MSUN 1.989e33 /* Solar mass */
#define YEAR 3.15576e7
#define MH 1.67352e-24 /* Mass of hydrogen atom. */
#define H 6.6260755e-27 /* Planck constant */
#define KB 1.38e-16 /*Boltzmann constant */
#define PI 3.141592654 /* PI - obviously. */
int npackets;
double lambda_min, lambda_max;
double freq_min, freq_max, freq_0;
double r0, r1;
double v1;
gsl_rng *rng; /* pointer for random number generator */
unsigned long int zseed; /* rnum generator seed */
struct spec
{
double lower_freq[MBINS];
double delta_freq[MBINS];
double flux[MBINS];
int nstart[MBINS];
} spectrum;
Makefile
CC = gcc
INCLUDE=/afs/mpa/common/pdsoft/include/
# use pg when you want to use gprof the profiler
#CFLAGS = -g -pg -Wall -I$(INCLUDE)
# Use this for large runs
CFLAGS = -O3 -Wall -I$(INCLUDE)
LIB=/afs/mpa/common/pdsoft/lib/
LDFLAGS= -L$(LIB) -lgsl -lgslcblas -lm
sl_objects = sl.o
sl: $(sl_objects)
$(CC) $(CFLAGS) $(sl_objects) $(LDFLAGS) -o sl