/* Code to calculate the spectrum of primordial magentic fields created by cosmic string loops in the early universe. If you use this code, please cite Diana Battefeld, Thorsten Battefeld, Daniel Wesley and Mark Wyman. "Magnetogenesis from cosmic string loops." arXiv:0708.2901 JCAP02:001 (2008) where more information about the physics encapsulated in the code can be found. The program tracks cohorts of loops that were created at similar times. It tracks the magnetic fields created on a variety of comoving scales. There are various dynamical processes that affect the properties of the loop cohorts which can be included. These dynamical processes can be also be switched off if desired. There are also two choices of string network models. The contribution from long strings is also calculated. Right now it compiles from the command line as gcc -o magnetiCS -lm magnetiCS.c To see a list of command-line options (which allow one to change network model or switch on/off various dynamical processes) do ./magnetiCS --help After the code finishes running, some diagnostic information is printed to standard error. The code also creates a file called "bfield_spectrum.txt" for output. The columns in this file are: column 1 j bin index column 2 u.B_loops[j].len_c comoving length for this bin (m) column 3 u.B_loops[j].a2B value of a^2 B (Gauss) column 4 u.B_loops[j].vol_c comoving volume covered by bin (m^3) column 5 u.B_longs[j].len_c comov. length for long string B (m) column 6 u.B_longs[j].a2B a^2 B for long strings (Gauss) You can also choose to output loop dynamics for an individual cohort. Use the option --follow-loops z1,z2,z3 to output information about the loop cohorts created at these three redshifts to "loop_dynamics.txt" This file contains the step number and redshift, followed by the physical length (m), transverse dimensionless velocity and rotational dimensionless velocity for each desired loop cohort. Version 1.0 */ /* DEFINES */ #define NUM_LOOP_COHORTS 1000 /* number of loop cohorts this needs to be a multiple of the num of steps in a, for the binning. If you adjust this variable, then adjust A_STEPS_PER_COHORT as well */ #define A_STEPS_PER_COHORT 1000 /* the number of steps in a whose loops are all binned together. This should satisfy NUM_LOOP_COHORTS * A_STEPS_PER_COHORT = NUM_A_STEPS */ #define NUM_A_STEPS 1000000 /* number of time steps to take. This must be a multiple of NUM_LOOP_COHORTS for the binning to work properly */ #define NUM_BFIELD_BINS 1000 #define NUM_A_EFOLDS 20 /* number of efoldings of a */ #define LN_A_STEP ((double)NUM_A_EFOLDS/(double)NUM_A_STEPS) /* log growth in a each setp */ /* Some cosmological parameters */ #define TODAY_H0 2.26854354e-18 /* in Hz, for h=.7 */ #define TODAY_OM 0.266 #define TODAY_OR (0.266 / 3454.0) /* use z_eq */ #define TODAY_OL (1.0-TODAY_OM-TODAY_OR) /* assume flat */ #define Z_DEC 1089 /* redshift at decoupling*/ #define Z_MR 3454 /* redshift at Matter-Radiation Equality*/ /* Physical constants */ #define EPIGO3 5.59035941e-10 /* 8\pi G/3 m^3/kg/s^2 */ #define M_PER_PC 3.08568025e16 /* meters per parsec */ #define SPEED_OF_LIGHT 2.99792458e8 /* m/s */ #define C2OG 1.34685326e27 /* c^2/G in kg/m */ #define S_PER_YEAR 3.15569260e7 /* seconds in a year */ #define S_PER_GYR 3.15569260e16 #define G_N 6.67300000e-11 /* Newton's G */ #define MP_E 0.000208793678 /* (M_proton / e) in Gauss - seconds */ /* String Network Constants */ #define C1RAD 0.21 /* model parameter c1 in the radiation era */ #define C2RAD 0.18 /* model parameter c2 in the radiation era */ #define C1MAT 0.2475 /* c1 in the matter (and dark energy) era(s) */ #define C2MAT 0.3675 /* c2 in the matter (and dark energy) era(s) */ #define P 0.28 /* string density parameter, global (for now) */ #define wigglyRAD 1.9 /* string wiggliness in RD */ #define wigglyMAT 1.5 /* string wiggliness in MD */ #define Gmu 2.e-7 /* string tension */ #define alph 1.e-2 /* new loop length as function of horizon size*/ /* minimum correlation length of magnetic fields, in meters, today (at z=0, a=1)*/ #define BFIELD_L_MIN 300 * M_PER_PC /* Let's say .3 kpc */ /* Now let's say we want to get between .3 kpc and 3 Gpc today. This is a factor of 1e7, whose natural log is 16.11809565095832. So the e-folding per bin would be */ #define BFIELD_EFOLD_PER_BIN (16.11809565095832/(double)NUM_BFIELD_BINS) /* INCLUDES */ #include #include #include #include #include #include /* For giving the user messages */ #define MSG(x) fprintf(stderr,x) #define MSG2(x,y) fprintf(stderr,x,y) #define MSG3(x,y,z) fprintf(stderr,x,y,z) FILE * so = stdout; /* The following options all have to do with global variables that are set by command-line options. */ /* String loop network mode */ #define OPTION_VOS 1 #define OPTION_OSM 2 /* String loop dynamics */ #define OPTION_LOOP_DYN_ON 1 #define OPTION_LOOP_DYN_OFF 0 #define LOOP_OUTPUT_FNAME "loop_dynamics.txt" #define BLOOP_OUTPUT_FNAME "bfield_spectrum.txt" #define OPT_NUM_HELP -2 #define OPT_NUM_NONE -1 /* Unsupported options */ #define OPT_NUM_OSM 0 #define OPT_NUM_VOS 1 #define OPT_NUM_LOOP_DYN_ON 2 #define OPT_NUM_LOOP_DYN_OFF 3 #define OPT_NUM_B_LOOPS_FILE 4 #define OPT_NUM_B_LONGS_FILE 5 #define OPT_NUM_COSMO_FILE 6 #define OPT_NUM_FOLLOW_LOOPS 7 #define OPT_NUM_FOLLOW_LOOPS_FILE 8 #define OPT_NUM_PRINT_STEPS 9 /* All possible command line options -- see man page for getopt_long() */ struct option opts[] = { { "help", no_argument, NULL, OPT_NUM_HELP }, /* Network model */ { "osm", no_argument, NULL, OPT_NUM_OSM }, { "vos", no_argument, NULL, OPT_NUM_VOS }, /* Loop dynamics */ { "loop-dyn", no_argument, NULL, OPT_NUM_LOOP_DYN_ON }, { "no-loop-dyn", no_argument, NULL, OPT_NUM_LOOP_DYN_OFF }, /* Output files for loop or long-string stuff. STDOUT will go to stdout */ { "B-loops-file", required_argument, NULL, OPT_NUM_B_LOOPS_FILE }, { "B-longs-file", required_argument, NULL, OPT_NUM_NONE }, /* Would you like some cosmology with that? */ { "cosmo-file", required_argument, NULL, OPT_NUM_NONE }, /* Follow individual cohorts */ { "follow-loops", required_argument, NULL, OPT_NUM_FOLLOW_LOOPS }, { "follow-loops-file", required_argument, NULL, OPT_NUM_FOLLOW_LOOPS_FILE }, /* How many points should be printed? (first and last are always printed) */ { "print-steps", required_argument, NULL, OPT_NUM_PRINT_STEPS }, /* End-of-array */ { NULL, 0, NULL, 0 } }; /* Help strings for each option. Keyed to the same index as in the option array above. So if you sandwich new options in, be sure to sandwich some help string in here too! (Even if it's just the empty string) */ char * opts_help[] = { "Help message", "Use the OSM string network model", "Use the VOS string network model (default)", "Turn on loop dynamics (default)", "Turn off loop dynamics", "Set the file for B-fields due to loops", "Set the file for B-fields due to long strings", "Set the file for output of cosmology info", "Follow loop cohorts created at z1,z2,z3", "Output file for loop cohort data", "Total number of steps that should be printed", NULL }; /* The number of special loop cohorts to follow, don't confuse with the total number of loop cohorts! */ #define MAX_COHORTS_TO_FOLLOW 20 typedef struct { int network_model; int loop_dynamics; int num_loop_follow; long loop_follow[MAX_COHORTS_TO_FOLLOW]; long print_block_size; char * loop_out_fname; FILE * loop_out; char * bloop_out_fname; FILE * bloop_out; } t_options; /* A global holding all the options! */ t_options g_opt; /* Nifty defines */ #define SQUARE(x) ((x)*(x)) /* Saves a pow call */ #define CUBE(x) ((x)*SQUARE(x)) /* ditto */ /* Program constants */ #define STEP_UNIVERSE_CONTINUE 1 #define STEP_UNIVERSE_FINISHED -1 /* Parameters */ const double PGROW= -0.20; /* parameters for Runge-Kutta integrator */ const double PSHRNK = -0.25; /* of standard VOS equations */ const double FCOR = 0.0666666666667; /* They have no physical meaning,*/ const double SAFETY = 0.9; /* but are taken from Numerical Recipes */ const double ERRCON = 6.0e-4; /* see www.nr.com for more */ const double hmax = 0.2; const double TINY = 1.0e-30; const double etamax = 20000; /* total amount of conformal time for universe to evolve. conformal time here is defined as deta = v dt / a */ const double starteta = 0.001; /* starting value for conformal time */ const double vi = 0.6; /* initial string RMS dimensionless velocity */ const double Hi = 1.0; /* initial Hubble constant */ const double era = 2.0; /* designates cosmological era: 2 for RD, 1.5 for MD */ const double PI = 3.14159265; /* the ratio of the circumference of circle to its diameter */ const double c1 = 0.21; /* one of the two free parameters in the VOS code for matching to simulations */ const double c2 = 0.18; /* the other */ /* TYPEDEFS */ /* NAME: t_loop_cohort Keeps track of a cohort of loops that were all created at nearly the same time. These loops will all have the same properties as the universe evolves, according to the one scale model. The lengths and so forth here are all quoted in physical (eg, not comoving units) as this is more convenient for describing the loop physics. To remind the programmers of this, a "p" is appended after the variable name. len_p Minimum length of loops in the cohort, m (physical units) d_len_p Spread of lengths in the cohort, m (physical units) u_t Translational velocity of loops, fraction of c u_r Rotational velocity of loops, fraction of c num Number density of loops in the cohort in comoving units. wig wiggliness of particular loop, determined by epoch of formation. */ typedef struct { double len_p; double d_len_p; double u_t; double u_r; double num_c; double wig; } t_loop_cohort; typedef struct { t_loop_cohort cohorts[NUM_LOOP_COHORTS]; double Gamma; double dL_dt; double Gamma_p; /* Describes rocket effect of GW emission on loops */ double Gamma_gr; /* Describes the torque due to GW emission */ double fr; /* The loop energy redshifting factor */ } t_loop_pop; /* NAME: t_bfield_bin Describes the strength of the magnetic field on a certain comoving scale. It is assumed that the magnetic fields, once created, will expand along with the expansion of the universe. Thus the units quoted here are comoving units. To remind the poor nitwitted programmers, the letter "c" is appended to the variable names. Furthermore, MHD says that B scales like (1+z)^2. Therefore, in order to avoid continuously correcting B for changes in the scale factor, this struct keeps track of B/(1+z)^2, or put another way the flux of B through a patch of constant comoving area. len_c The length scale of the magnetic field (comoving units) d_len_c Width of the bin (comoving units) a2B Magnetic field flux vol_c The volume filled by the magnetic field, as a fraction of the current horizon, for weighting purposes */ typedef struct { double len_c; double d_len_c; double a2B; double vol_c; } t_bfield_bin; /* NAME: t_univ_state Keeps track of the state of the universe t The time since the Big Bang (seconds) d_t Time step to the next step in a. a Scale factor (1=today) z Redshift (0=today) rho_r Energy density in radiation (kg/m^3) rho_m Energy density in matter (baryonic+cdm) rho_l Energy density in cosmological constant rho Total density (kg/m^3) H Hubble parameter (in Hz) ph Particle horizon (m) physical units d_ph Change in the particle horizon from the next step in a (m) rho_ls Energy density in long strings (not actually kept here) vol_loop Volume of space that has been swept by loops dvol New Volume swept in last time step c1 VOS parameter c1 c2 VOS parameter c2 */ typedef struct { double t; double d_t; double a; double z; double rho_r; double rho_m; double rho_l; double rho; double H; double ph; double d_ph; double vol_loop; double dvol; double c1; double c2; double wiggliness; } t_univ_state; /* NAME: t_long_string_pop Keeps track of quantities related to the long string population rho Energy density in long strings Gmu_c2 G mu / c^2 (dimensionless) mu in kg/m d_rho_dt Rate of change of rho with cosmic time A Appears in energy density formula for long strings mean_u2 Appears in loop creation rate formula Lstraight_c comoving string "straightness" length in hubble units; multiply times a and the l_ph to get meters Lcorr_c =1/comoving number density, the average comoving distance between strings, in hubble units Lstraight_p same as above, but Lcorr_p put into physical mks units */ typedef struct { /* These shouldn't change over time */ double Gmu_c2; double mu; double alpha; /* These change over time */ double rho; double d_rho_dt; double Lstraight_c; double Lcorr_c; double Lstraight_p; double Lcorr_p; /* The following parameters have to do with the scaling assumption that the energy density in long strings goes like \rho = A \mu c^2 / ph^2 with A a constant that changes depending on whether we are in MD or RD. We also have a mean-squared rapidity u (pure) which also changed in a way that depends on the era. */ double A; double A_MD; double A_RD; double mean_u2; double mean_u2_RD; double mean_u2_MD; } t_long_string_pop; /* NAME: t_universe Keeps everything together. loop The loop cohorts longs The long string population B_loops Magnetic field bins for loops B_longs Magnetic field bins for long strings univ Current state of the universe step Which step we're in of NUM_A_STEPS */ typedef struct { t_loop_pop loops; t_long_string_pop longs; t_bfield_bin B_loops[NUM_BFIELD_BINS]; t_bfield_bin B_longs[NUM_BFIELD_BINS]; t_univ_state state; long step; } t_universe; /* DECLARATIONS */ /* Initialization functions */ void init_options( t_options * o , int argc, char * argv [] ); void init_options_open_file( FILE ** f, char * name ); void init_options_help( void ); void init_options_help_one( int j ); void init_universe( t_universe * u ); void init_loops( t_universe * u ); void init_B( t_universe * u ); void init_univ_state( t_universe * u ); void init_long_string_pop( t_universe * u ); /* Stepper functions */ int step_universe( t_universe * u ); void step_long_string_pop( t_universe * u_new, t_universe * u_old ); void step_loops( t_universe * u_new, t_universe * u_old ); void step_univ_state( t_universe * u_new, t_universe * u_old ); void step_B( t_universe * u_new, t_universe * u_old ); /* VOS equations solvers */ double f(double c1, double vp, double Ha); double g(double c2, double v, double Ha, double ell); void rk4( double * y, double * dy, double * yout); double hfunc(double c2, double v, double Ha, double ell, double sL); /* Printing data structures */ void print_universe( FILE * f, t_universe * u ); void print_long_string_pop( FILE * f, t_long_string_pop *l ); void print_loop_pop( FILE * f, t_loop_pop * lp ); void print_loop_cohort( FILE * f,t_loop_cohort * lc ); void print_univ_state( FILE * f, t_univ_state * us ); void univ_state_at_a( t_univ_state * us, double a ); void univ_state_at_j( t_univ_state * us, long j ); void long_string_pop( t_long_string_pop * ls, t_univ_state * us ); void long_string_pop_scaling_ansatz( t_long_string_pop *ls, \ t_univ_state * us ); void loops_add( t_universe * u ); int main( int argv, char * argc[] ) { long j; double vol; t_universe u; long count = 0; long cohort; /* Initialise the options */ init_options( &g_opt, argv, argc ); /* Fiat lux! */ init_universe( &u ); /* Open the loop data file if not opened yet */ if ( g_opt.num_loop_follow > 0 ) init_options_open_file( &(g_opt.loop_out), LOOP_OUTPUT_FNAME ); init_options_open_file( &(g_opt.bloop_out), g_opt.bloop_out_fname ); /* Continue on as long as we're supposed to. */ while ( step_universe( &u ) == STEP_UNIVERSE_CONTINUE ) { /* Only print when we're supposed to! */ if ( (count % g_opt.print_block_size == 0) || \ (count == NUM_A_STEPS-1) ) { /* First: are we tracking any loops? */ if ( g_opt.num_loop_follow > 0 ) { fprintf( g_opt.loop_out, "%10ld %+1.5e ", count, u.state.z ); for ( j=0; j< g_opt.num_loop_follow; j++ ) { cohort = g_opt.loop_follow[j]; fprintf( g_opt.loop_out, "%+1.5e %+1.5e %+1.5e ", u.loops.cohorts[cohort].len_p, u.loops.cohorts[cohort].u_t, u.loops.cohorts[cohort].u_r ); } fprintf( g_opt.loop_out, "\n" ); } } count++; } fprintf( stderr, "Current age: %e s = %e Gyr\n", \ u.state.t, u.state.t / S_PER_GYR ); print_universe( stderr, &u ); /* there are about 2.3e78 m^3 per horizon */ vol = 0; for (j=0; j 0) { fprintf( g_opt.bloop_out, "%10d %+1.5e %+1.5e %+1.5e %+1.5e %+1.5e\n", j, u.B_loops[j].len_c, u.B_loops[j].a2B, u.B_loops[j].vol_c, u.B_longs[j].len_c, u.B_longs[j].a2B ); vol += u.B_loops[j].vol_c; } } /* Oddly, this'll close STDOUT. */ if ( g_opt.loop_out != NULL ) fclose(g_opt.loop_out); if ( g_opt.bloop_out != NULL ) fclose(g_opt.bloop_out); return 0; } /* NAME: init_options ARGS: argc Number of command-line arguments argv Array of strings for command-line tokens o Pointer to options struct. RETS: void This function takes care of reading command line options into the code. */ void init_options( t_options * o , int argc, char * argv [] ) { int j, opt_num; extern char * optarg; char * tok, c; double cohort_a, cohort_z; long cohort_j; /* Give the options array the proper number in which to return vals */ for ( j=0; opts[j].name != NULL; j++ ) opts[j].flag = &opt_num; /* Set default options */ o->network_model = OPTION_VOS; o->loop_dynamics = OPTION_LOOP_DYN_ON; o->num_loop_follow = 0; o->print_block_size = NUM_A_STEPS/1000; o->loop_out_fname = LOOP_OUTPUT_FNAME; o->loop_out = NULL; o->bloop_out_fname = BLOOP_OUTPUT_FNAME; o->bloop_out = NULL; /* Mode where it just translates to short options */ c = getopt_long( argc, argv, "", opts, NULL); while ( c != -1 ) { /* Capture malformed argument lists. */ if ( (c=='?') || (c==':') ) { MSG( "I don't understand the syntax.\n" ); exit(-1); } switch (opt_num) { case OPT_NUM_HELP: init_options_help(); break; case OPT_NUM_NONE: MSG( "options: sorry, one option you sent is currently unsupported.\n" ); break; case OPT_NUM_VOS: MSG( "options: setting to VOS model\n" ); o->network_model = OPTION_VOS; break; case OPT_NUM_OSM: MSG( "options: setting to OSM model\n" ); o->network_model = OPTION_OSM; break; case OPT_NUM_LOOP_DYN_ON: MSG( "options: setting loop dynamics ON\n" ); o->loop_dynamics = OPTION_LOOP_DYN_ON; break; case OPT_NUM_LOOP_DYN_OFF: MSG( "options: setting loop dynamics OFF\n" ); o->loop_dynamics = OPTION_LOOP_DYN_OFF; break; case OPT_NUM_FOLLOW_LOOPS: MSG( "options: following some loops\n" ); break; case OPT_NUM_FOLLOW_LOOPS_FILE: o->loop_out_fname = optarg; MSG2( "options: loop data will go into \"%s\"\n", o->loop_out_fname); init_options_open_file( &(o->loop_out), optarg ); break; case OPT_NUM_B_LOOPS_FILE: o->bloop_out_fname = optarg; MSG2( "loop B-field data will go into \"%s\"\n", o->bloop_out_fname ); init_options_open_file( &(o->bloop_out), optarg ); break; case OPT_NUM_PRINT_STEPS: o->print_block_size = NUM_A_STEPS/atol( optarg ); MSG2( "will print every %d steps\n", o->print_block_size ); break; default: MSG( "wtf? unknown option!\n" ); break; } /* Some options require custom processing: If we have the loop redshift stuff, then tokenize! The user has requested that we follow some loops. We should oblige them by getting the requested loop creation redshifts, converting them to cohorts, saving this information to the options struct. */ if ( opt_num == OPT_NUM_FOLLOW_LOOPS ) { /* Tokenize the string of loop redshifts. */ tok = strtok( optarg, "," ); /* Read through the tokens. Each should be a float. */ while ( tok != NULL && \ (o->num_loop_followloop_follow[o->num_loop_follow] = cohort_j; MSG3( "options: ...following loops made at z=%lf cohort %ld\n", cohort_z, o->loop_follow[o->num_loop_follow] ); /* Increment the number of loop cohorts we're following. */ (o->num_loop_follow)++; /* Get next token */ tok = strtok( NULL, "," ); } } /* Get the next option (if any) */ c = getopt_long( argc, argv, "", opts, NULL); } } /* NAME: init_option_open_file ARGS: f Handle to a FILE struct name Name of file to open RETS: void Helper to open a file. Has a few behaviors designed to make this a useful replacement for a few lines of code that we have: 1. If *f is not NULL, it assumes the file has already been opened and does nothing (lets you call it with a default value which one uses if the file isn't already opened). 2. If name == STDOUT or stdout it sets *f to be stdout. */ void init_options_open_file( FILE ** f, char * name ) { assert( name != NULL ); if ( *f != NULL ) { MSG2( "file opener: declined to open \"%s\" as *f != NULL\n", name ); return; } /* Is it STDOUT? */ if ( (strcmp( "STDOUT", name ) == 0) || \ (strcmp( "stdout", name ) == 0) ) { MSG( "file opener: using stdout\n" ); *f = stdout; } else { MSG2( "file opener: opening \"%s\"\n", name ); *f = fopen( name, "w" ); assert( *f != NULL ); } } void init_options_help( void ) { int j; /* Supported options */ MSG( "\nSupported options:\n[arg] denotes a required argument.\n\n" ); for ( j=0; opts[j].name != NULL; j++ ) if ( opts[j].val != OPT_NUM_NONE ) init_options_help_one(j); /* Un-Supported options */ MSG( "\nUnsupported options:\n\n" ); for ( j=0; opts[j].name != NULL; j++ ) if ( opts[j].val == OPT_NUM_NONE ) init_options_help_one(j); /* Extra info. */ MSG( "\nThe filename STDOUT or stdout goes to the screen.\n\n"); exit(-1); } void init_options_help_one( int j ) { MSG2( "--%-20s ", opts[j].name ); if ( opts[j].has_arg != no_argument ) MSG ( "[arg] " ); else MSG ( " " ); MSG2( " %s\n", opts_help[j] ); } /* NAME: init_universe ARGS: u Pointer to Universe struct RETS: void This function initializes the universe in some initial state. The idea is that, since we know that log(a)=0 today, the time steps will be determined by taking equally spaced intervals in d log(a) until we reach the present, with the number of such steps determined by NUM_A_STEPS. */ void init_universe( t_universe * u ) { /* initialize the state of the universe */ init_univ_state( u ); /* initialize long string population */ init_long_string_pop( u ); /* initialize loop cohorts */ init_loops( u ); /* initialize bfield stuff */ init_B( u ); /* First step */ u->step = 0; } void init_loops( t_universe * u ) { long j; /* Now add lots of empty bins */ for (j=0; jloops.cohorts[j].len_p = 0.0; u->loops.cohorts[j].d_len_p = 0.0; u->loops.cohorts[j].num_c = 0.0; u->loops.cohorts[j].u_t = 0.1; u->loops.cohorts[j].u_r = 0.4; } u->loops.Gamma = 50.0; /* Because of the loop physics, the loss of length per time is the same for all the loops \dot{L} = Gamma * (G \mu / c^2) * c */ u->loops.dL_dt = u->loops.Gamma * u->longs.Gmu_c2 * SPEED_OF_LIGHT; u->loops.Gamma_p = 10.; u->loops.Gamma_gr = 5.; u->loops.fr = 0.71; } /* NAME: init_B ARGS: u: the Universe RETS: void Initialize the magnetic field arrays for the long strings and the loops. */ void init_B( t_universe * u ) { long j; /* Now add lots of empty bins */ for (j=0; jB_loops[j].len_c = 0.0; u->B_loops[j].d_len_c = 0.0; u->B_loops[j].a2B = 0.0; u->B_loops[j].vol_c = 0.0; u->B_loops[j].d_len_c = BFIELD_L_MIN*exp((double)j*BFIELD_EFOLD_PER_BIN); } for (j=0; jB_longs[j].len_c = 0.0; u->B_longs[j].d_len_c = 0.0; u->B_longs[j].a2B = 0.0; u->B_longs[j].vol_c = 0.0; u->B_longs[j].d_len_c = BFIELD_L_MIN*exp((double)j*BFIELD_EFOLD_PER_BIN); } } void init_univ_state( t_universe * u ) { u->state.t = 0.0; /* a kludge */ u->state.ph = 1. ; /* also, a kludge */ u->state.vol_loop = 0.0; /* not a kludge*/ u->state.dvol = 0.0; u->state.a = exp( -1.0 * NUM_A_EFOLDS ); univ_state_at_a( &(u->state), u->state.a ); } void init_long_string_pop( t_universe * u ) { u->longs.Gmu_c2 = Gmu; u->longs.mu = (u->longs.Gmu_c2) * C2OG; u->longs.alpha = alph; if( g_opt.network_model == OPTION_OSM ) { u->longs.Lstraight_p = 1.0; u->longs.Lcorr_p = 1.0 / sqrt(52.0); } if( g_opt.network_model == OPTION_VOS) { fprintf( stderr, "We're using VOS!\n" ); u->longs.Lstraight_c = 0.16/(u->state.H*u->state.a); u->longs.Lcorr_c = 0.0005/(u->state.H*u->state.a); u->longs.Lstraight_p = 0.16/u->state.H; u->longs.Lcorr_p = 0.0005/u->state.H; u->longs.mean_u2 = 0.43; } u->longs.A_RD = 52.0; u->longs.A_MD = 31.0; u->longs.mean_u2_RD = 0.43; u->longs.mean_u2_MD = 0.37; if(g_opt.network_model == OPTION_OSM) { long_string_pop( &(u->longs), &(u->state) ); } } /* NAME: long_string_pop ARGS: ls Pointer to a long string population struct us Pointer to a universe state struct RETS: void Given the state of the universe, this function determines features of the long string population, subject to the scaling assumption. */ void long_string_pop( t_long_string_pop * ls, t_univ_state * us ) { int temp; t_univ_state us_2; /* spare universes and long string pops */ t_long_string_pop ls_2; double new_a; /* Fix properties of the long string population according to scaling */ long_string_pop_scaling_ansatz( ls, us ); /* Now we need to calculate the change in rho with time. We can do this because d ln(rho) / dt = -2 d ln(ph) /dt + (terms involving change in weights). If we ignore the latter ones, since they're going to be mostly small, are only important near M/R transition, we will get accumulated errors in teh total loops produced. So let's factor them in as well. To do this all properly, copy the structs over and then increment the time by one to allow the routines controlling the universe to get everything right. Then call the helper function to calculate the new rho, then finally just subtract to get the derivative! Make copies of the long string population and universe in its current state */ memcpy( &us_2, us, sizeof(t_univ_state) ); memcpy( &ls_2, ls, sizeof(t_long_string_pop) ); /* Increment cosmic time by a single step in a */ new_a = us->a * exp( LN_A_STEP ); univ_state_at_a( &us_2, new_a); us_2.ph += us_2.d_ph; /* Calculate rho in the new universe one time step in the future */ long_string_pop_scaling_ansatz( &ls_2, &us_2 ); ls->d_rho_dt = (ls_2.rho - ls->rho)/(us->d_t); ls->Lstraight_p = us->ph; ls->Lcorr_p = us->ph / sqrt(ls->A); } void long_string_pop_scaling_ansatz( t_long_string_pop *ls, \ t_univ_state * us ) { double r_weight, m_weight; /* Radiation or matter domination? To prevent numerical issues, we smoothly interpolate between the two regimes, by weighting A and mean_u2 according to the densities in each component. This avoids divergent loop creation rates as the time step goes to zero. After MD, we just use the MD values -- this holds even when we are deep in the dark-energy dominated regime. The long string energy density formula is taken from Caldwell, Allen PRD _45_ 3447 (3.1). Right now there is a dimensionless parameter "A" which appears in this formula, and we smoothly interpolate between the RD and MD values. */ r_weight = us->rho_r / us->rho; m_weight = 1.0 - r_weight; ls->A = (r_weight * (ls->A_RD)) + (m_weight * (ls->A_MD)); ls->mean_u2 = r_weight * (ls->mean_u2_RD) + m_weight * (ls->mean_u2_MD); /* Apply scaling ansatz */ /* rho is a MASS PER VOLUME, not an ENERGY PER VOLUME so that Friedmann works ok */ ls->rho = (ls->A) * pow(SPEED_OF_LIGHT,2) * (ls->mu) *pow(us->ph,-2.0); } /* NAME: loops_add ARGS: RETS: void Applies the OSM assumptions to compute the number of loops that are created. These are then added to the loop cohort arrays in the appropriate place. */ void loops_add( t_universe * u ) { double n, len, w, co_w; long j; /* Equation (3.7) in Caldwell and Allen, PRD 45 (1992) 3447 */ n = -1.0 * pow(u->state.a, 3.0 ) / ( u->longs.mu * u->longs.alpha * \ u->state.ph * SPEED_OF_LIGHT * SPEED_OF_LIGHT ); n *= u->longs.d_rho_dt + 2.0 * u->state.H * u->longs.rho * \ (1.0 + u->longs.mean_u2 ); n *= u->state.d_t; /* The length of the newly created loops */ len = u->longs.alpha * u->state.ph * u->loops.fr; /* n is now the number density of strings created during this time step dt. Next we need to put them in the proper bin. */ j = u->step / A_STEPS_PER_COHORT; /* bounds checking */ if ( (j<0) || (j>(NUM_LOOP_COHORTS-1)) ) { fprintf( stderr, " uh oh! \n" ); return; } /* add the loops to the loop cohorts. w is a weight, so that for example if we are adding 10 loops with average velocity a to a cohort which already has 90 loops with average velocity b, then new average velocity will be (10*a + 90*b)/100. */ w = n/(n + u->loops.cohorts[j].num_c); co_w = 1.0 - w; if(n<0) {n = 0; w=0.5; co_w=0.5;} /* weight physical properties by the number of loops */ u->loops.cohorts[j].num_c += n; u->loops.cohorts[j].len_p = (w*len + co_w * u->loops.cohorts[j].len_p ); u->loops.cohorts[j].wig = (w*u->state.wiggliness + co_w * u->loops.cohorts[j].wig ); } /* NAME: univ_step ARGS: RETS: STEP_UNIVERSE_FINISHED when it's done STEP_UNIVERSE_CONTINUE when not. Iterates the universe a single step in time. */ int step_universe( t_universe * u ) { t_universe u_new; if ( u->step >= (NUM_A_STEPS-1) ) return STEP_UNIVERSE_FINISHED; memcpy( &u_new, u, sizeof(t_universe) ); step_univ_state( &u_new, u ); step_long_string_pop( &u_new, u ); step_loops( &u_new, u ); step_B( &u_new, u); memcpy( u, &u_new, sizeof(t_universe) ); (u->step)++; return STEP_UNIVERSE_CONTINUE; } void step_univ_state( t_universe * u_new, t_universe * u_old ) { /* increment the time and particle horizon */ u_new->state.t = u_old->state.t + u_old->state.d_t; u_new->state.ph = u_old->state.ph + u_old->state.d_ph; u_new->state.vol_loop = u_old->state.vol_loop + u_new->state.dvol; /* everything else is taken care of here */ univ_state_at_j( &(u_new->state), u_old->step + 1 ); } void step_loops( t_universe * u_new, t_universe * u_old ) { long j; long current; double sum=0, num=0; double dt,len; double stringpot,lambda,Gu,Gu_c2; current = (u_new->step)/A_STEPS_PER_COHORT; dt = u_old->state.d_t; Gu = SQUARE(SPEED_OF_LIGHT)*u_old->longs.Gmu_c2; Gu_c2 = u_old->longs.Gmu_c2; /* new loops */ loops_add( u_new ); /* The newly created ones lead to a subtle bug, since they have zero physical length. */ u_new->state.dvol = 0; for (j=0; jloops.cohorts[j].num_c > 0.0) && (u_old->loops.cohorts[j].len_p > 0.0) && (u_new->loops.cohorts[j].len_p > 0.0)) { /* Take account of gravitational radiation */ /* Shrink the loops that are already created */ num = u_new->loops.cohorts[j].num_c; u_new->state.dvol += num*pow(u_new->state.a,-3.0)*pow(u_new->loops.cohorts[j].len_p,2)*u_new->loops.cohorts[j].u_t*SPEED_OF_LIGHT*u_old->state.d_t/(4*PI); sum += num; u_new->loops.cohorts[j].len_p -= u_new->loops.dL_dt * \ u_old->state.d_t; /* This takes care of the loops that are already gone */ if (u_new->loops.cohorts[j].len_p < 0.0) { u_new->loops.cohorts[j].num_c = 0.0; u_new->loops.cohorts[j].len_p = 0.0; u_new->loops.cohorts[j].u_t = 0.0; u_new->loops.cohorts[j].u_r = 0.0; } /* Physical length of the loops after all the mods */ len = u_new->loops.cohorts[j].len_p; stringpot = (u_new->loops.cohorts[j].wig - 1)*u_new->longs.Gmu_c2; lambda = stringpot * SQUARE(SPEED_OF_LIGHT) / G_N; /* Equations governing the loop velocities and such */ if (( g_opt.loop_dynamics == OPTION_LOOP_DYN_ON ) && (len>0.)) { /* Do the change to translational velocity */ double v_t; double rho; double t_star; double H_friction; double dyn_friction; double rocket; /* These have to use the NEW values since loops_add may do something. */ v_t = u_new->loops.cohorts[j].u_t * SPEED_OF_LIGHT; rho = u_old->state.rho; t_star = CUBE(v_t)/( 4.*PI*SQUARE(G_N)*len*lambda*rho ); H_friction = -1. * u_old->state.H * v_t; dyn_friction = (-1. * v_t / t_star) * ( log(1.5e4/ u_old->longs.alpha) ); rocket = Gu * u_old->loops.Gamma_p / (len); u_new->loops.cohorts[j].u_t = u_old->loops.cohorts[j].u_t \ + ((H_friction + dyn_friction + rocket)/SPEED_OF_LIGHT)*dt; /* Do the change in rotational velocity */ double torque_gw; double change_in_I; /* This is not actually the torque, but is due to the torque */ torque_gw = -1. * (SQUARE(Gu)/G_N) * u_old->loops.Gamma_gr /\ (lambda*len*SPEED_OF_LIGHT); change_in_I = 2. * u_old->loops.Gamma * Gu_c2 * \ u_old->loops.cohorts[j].u_r*SPEED_OF_LIGHT / len; u_new->loops.cohorts[j].u_r = u_old->loops.cohorts[j].u_r \ + (torque_gw + change_in_I)*dt; if ( u_new->loops.cohorts[j].u_r < 0.) u_new->loops.cohorts[j].u_r = 0.; if ( u_new->loops.cohorts[j].u_t < 0.) u_new->loops.cohorts[j].u_t = 0.; if ( u_new->loops.cohorts[j].u_r > 0.7071 ) u_new->loops.cohorts[j].u_r = 0.7071; if ( u_new->loops.cohorts[j].u_t > 0.7071 ) u_new->loops.cohorts[j].u_t = 0.7071; } } } } /* NAME: step_B ARGS: u_new Pointer to the new universe struct with all elements updated u_old Pointer to the old (current) universe struct RETS: void The stepper function for magnetic field generation. Probably one shouldn't have u_new and u_old pointing at the same struct. */ void step_B( t_universe * u_new, t_universe * u_old ) { long j; long current; int len_step; double sum=0, num=0, numold; double bmag, length, volume; double stringpot; /* G (T - mu) / c^2 */ /* the minimum length, redshifted to the current epoch */ double minLnow, binLnow; double weight=0; /* weighting factor */ double string_u, l_corr; double z, dvol_c; minLnow = BFIELD_L_MIN/(1+u_new->state.z); /* stringpot is G\lambda/c^2 */ stringpot = (u_new->state.wiggliness - 1)*u_new->longs.Gmu_c2; current = (u_new->step)/A_STEPS_PER_COHORT; z = u_new->state.z; /* Calculate B-fields and add them to the appropriate length bin */ /* Check to see that the universe is still ionized */ if((u_new->state.z>1000)) { /* Now for the contribution from the loops */ for (j=0; jloops.cohorts[j].num_c > 0.0) && \ (u_new->loops.cohorts[j].len_p > minLnow) && \ (u_new->loops.cohorts[j].u_t < u_new->loops.cohorts[j].u_r)) { numold = volume; num = u_new->loops.cohorts[j].num_c; /* string pot is G \lambda / c^2, for this loop cohort */ stringpot = (u_new->loops.cohorts[j].wig - 1)*u_new->longs.Gmu_c2; /* Compute the contribution to volume swept by this cohort */ volume = num * \ pow(u_new->state.a,-3.0) * \ pow(u_new->loops.cohorts[j].len_p,2) * \ u_new->loops.cohorts[j].u_t*SPEED_OF_LIGHT*u_old->state.d_t; /* Strength of the magnetic field */ bmag = SPEED_OF_LIGHT * (2./7.)*MP_E*4.*PI*PI* \ SQUARE(stringpot)/\ (u_new->loops.cohorts[j].len_p * u_new->loops.cohorts[j].u_r * \ SQUARE(u_new->loops.cohorts[j].u_t) ); /* Now multiply by a^2 in order to get the constant flux */ bmag = bmag / SQUARE( 1. + z ); length = u_new->loops.cohorts[j].len_p; /* len_step is the index in the bfield binning array. */ len_step = (int)( log(length/minLnow)/(BFIELD_EFOLD_PER_BIN*0.9) ); /* Put all the lengths that are too big in the last bin */ if(len_step>=NUM_BFIELD_BINS) len_step = NUM_BFIELD_BINS-1; /* Of course, what if len_step is too small? looking again, it seems that this case is alreay excluded by the if statement above, but why not just do it by symmetry */ if ( len_step < 0) len_step = 0; weight = volume + u_new->B_loops[len_step].vol_c; u_new->B_loops[len_step].a2B = (volume*bmag + u_new->B_loops[len_step].vol_c*u_new->B_loops[len_step].a2B)/weight; u_new->B_loops[len_step].vol_c = u_new->B_loops[len_step].vol_c + volume; u_new->B_loops[len_step].len_c = (volume*u_new->B_loops[len_step].len_c + u_new->B_loops[len_step].vol_c*(u_new->loops.cohorts[j].len_p/u_new->state.a))/weight; } } /* end stepping over loop cohorts */ /* Now for the contribution from the long strings References are Davis, Dimopoulos PRD 72, 043517 (2005) but we use the formulae from our paper here. First of all, the correlation length here is the inter-string distance We could use D+D (36) for this. However, from the arguments in D+D and the osm it seems better to use something like the particle horizon divided by N, where N is the (constant) number of strings per particle horizon expected by the one-scale model. */ l_corr = u_new->longs.Lcorr_p; /* Physical corr length */ /* long strings have longer intrinsic length scales, so a different binning makes more sense: */ minLnow *= 1; binLnow *= 50; if((l_corr>minLnow)) { /* Take R_s to be the correlation length */ string_u = sqrt(u_new->longs.mean_u2); /* stringpot here is the long string wiggliness, determined by epoch */ stringpot = (u_new->state.wiggliness - 1)*u_new->longs.Gmu_c2; bmag = SPEED_OF_LIGHT * (2.)*MP_E * 4 * PI*PI * SQUARE(stringpot) /\ (2.*l_corr * CUBE(string_u)); /* Modify to keep magnetic flux constant */ bmag = bmag / SQUARE( 1. + z ); /* Now put this contribution in a bin. One issue is the weighting factor: how should we weight the long string encounters? The length scale is determined by the inter-string distance only, so we should average over the b-fields produced by strings during the time when the correlation length is within the bin edges. Let's average over time. I will mimic this by using the vol_c member of the B-field bin struct. */ dvol_c = u_old->state.d_t * string_u* SPEED_OF_LIGHT; /* Choose which bin to put this in */ j = log(l_corr/minLnow)/BFIELD_EFOLD_PER_BIN; if(j>=NUM_BFIELD_BINS) j = NUM_BFIELD_BINS -1; /* Do a weighted addition of the field to the bins. */ weight = (dvol_c + u_new->B_loops[j].vol_c ); u_new->B_longs[j].a2B = (dvol_c * bmag + u_new->B_loops[j].vol_c*u_new->B_longs[j].a2B)/weight; /* Set the comoving length */ u_new->B_longs[j].len_c = \ (dvol_c * (1. + z ) * l_corr + u_new->B_loops[j].vol_c* u_new->B_longs[j].len_c)/weight; u_new->B_longs[j].vol_c += dvol_c; } } } void step_long_string_pop( t_universe * u_new, t_universe * u_old ) { double y[7], dy[6], yout[6]; double gamma, oldrho; if(g_opt.network_model == OPTION_OSM) { /* No integration */ long_string_pop( &(u_new->longs), &(u_new->state) ); } else { y[0] = u_new->state.c1; y[1] = u_new->state.c2; y[2] = u_new->state.H; y[3] = sqrt(u_new->longs.mean_u2); y[4] = u_new->longs.Lstraight_c; y[5] = u_new->longs.Lcorr_c; y[6] = u_new->state.a; rk4(y,dy,yout); u_new->longs.mean_u2 = SQUARE(yout[3]); u_new->longs.Lstraight_c = yout[4]; u_new->longs.Lcorr_c = yout[5]; u_new->longs.Lstraight_p = u_new->longs.Lstraight_c*u_new->state.a*u_new->state.H*u_new->state.ph; u_new->longs.Lcorr_p = u_new->longs.Lcorr_c*u_new->state.a*u_new->state.H*u_new->state.ph; gamma = 1/sqrt(1-u_new->longs.mean_u2); u_new->longs.rho = SQUARE(SPEED_OF_LIGHT)*u_new->longs.mu/SQUARE(u_new->longs.Lcorr_p); u_new->longs.d_rho_dt = (u_new->longs.rho - u_old->longs.rho)/u_new->state.d_t; } } /* NAME: univ_state_at_a ARGS: univ_state Pointer to a universe state struct a The desired scale factor RETS: void This function solves the Friedmann equation and associated things. For the value of "a" passed, this function computes the various densities and so on. It solves the Friedmann equation by also giving a value for the Hubble parameter. It also sets the two differential quantities dt and dph, which are the amount added to the time and particle horizon (respectively) between now and the next step in a. NOTE: the function does not compute t and ph, as these are integrated quantities! It does however compute d_t and d_ph. It also automatically computes the redshift (as a convenience). */ void univ_state_at_a( t_univ_state * us, double a ) { double scaled_or, /* scaled variables */ scaled_om, scaled_ol; double scaledc1, scaledc2, wiggly; /* These are today's Omega parameters, scaled by the growth of the relevant energy with scale factor. */ us->a = a; scaled_or = TODAY_OR * pow(a,-4.0); scaled_om = TODAY_OM * pow(a,-3.0); scaled_ol = TODAY_OL; scaledc1 = (C1RAD*scaled_or + C1MAT*scaled_om)/(scaled_or + scaled_om); scaledc2 = (C2RAD*scaled_or + C2MAT*scaled_om)/(scaled_or + scaled_om); wiggly = (wigglyRAD*scaled_or + wigglyMAT*scaled_om)/(scaled_or + scaled_om); us->z = 1.0 / us->a; us->H = TODAY_H0 * sqrt( scaled_om + scaled_or + scaled_ol ); us->rho_r = (TODAY_H0*TODAY_H0/EPIGO3)*scaled_or; us->rho_m = (TODAY_H0*TODAY_H0/EPIGO3)*scaled_om; us->rho_l = (TODAY_H0*TODAY_H0/EPIGO3)*scaled_ol; us->rho = us->rho_r + us->rho_m + us->rho_l; us->c1 = scaledc1; us->c2 = scaledc2; us->wiggliness = wiggly; /* Compute the time step till the next step in a. Since H = d ln(a) / dt, we know that dt = d ln(a) / H. */ us->d_t = LN_A_STEP / us->H; /* The change in the particle horizon comes from the expansion of the current particle horizon, plus some additional from the further distance traveleed by the light ray in the time dt. Below, we will assume that LN_A_STEP is very small. This is to avoid having to do e^x - 1 which will lead to round-off errors */ assert( LN_A_STEP*LN_A_STEP*LN_A_STEP < 1e-9 ); us->d_ph = us->ph * LN_A_STEP * ( 1.0 + 0.5 * LN_A_STEP ) + \ SPEED_OF_LIGHT * us->d_t; } /* NAME: univ_state_at_j ARGS: us Pointer to a universe state struct j which step Given a step 0 < j < NUM_A_STEPS-1, this gives the state of the universe at this step. */ void univ_state_at_j( t_univ_state * us, long j ) { us->a = exp( -1.0*NUM_A_EFOLDS*((double)(NUM_A_STEPS-j-1))/((double)(NUM_A_STEPS-1)) ); univ_state_at_a( us, us->a ); } void print_universe( FILE * f, t_universe * u ) { fprintf( f, "universe.step %ld\n", u->step ); print_univ_state(f, &(u->state) ); print_long_string_pop(f, &(u->longs) ); print_loop_pop(f, &(u->loops) ); } void print_loop_cohort( FILE * f,t_loop_cohort * lc ) { fprintf( f, "\t\tlen_p %e\n", lc->len_p ); fprintf( f, "\t\td_len_p %e\n", lc->d_len_p ); fprintf( f, "\t\tu_t %e\n", lc->u_t ); fprintf( f, "\t\tu_r %e\n", lc->u_r ); fprintf( f, "\t\tnum_c %e\n", lc->num_c ); } void print_loop_pop( FILE * f,t_loop_pop *lp ) { fprintf( f, "\n\tloop_pop:\n" ); fprintf( f, "\tGamma %e\n", lp->Gamma ); fprintf( f, "\tdL_dt %e\n", lp->dL_dt ); } void print_univ_state( FILE * f,t_univ_state *us ) { fprintf( f, "\n\tuniv_state:\n" ); fprintf( f, "\tt %e = %lf Gyr\n", us->t, us->t / S_PER_GYR ); fprintf( f, "\td_t %e = %lf yr\n", us->d_t, us->d_t / S_PER_YEAR ); fprintf( f, "\ta %e\n", us->a ); fprintf( f, "\tz %e\n", us->z ); fprintf( f, "\trho_r %e\n", us->rho_r ); fprintf( f, "\trho_m %e\n", us->rho_m ); fprintf( f, "\trho_l %e\n", us->rho_l ); fprintf( f, "\trho %e\n", us->rho ); fprintf( f, "\tH %e\n", us->H ); fprintf( f, "\tph %e = %e Glyr\n", us->ph, us->ph / (S_PER_YEAR * SPEED_OF_LIGHT * 1e9) ); fprintf( f, "\td_ph %e\n", us->d_ph ); fprintf( f, "\tvol_loop %e\n", us->vol_loop); } void print_long_string_pop( FILE * f, t_long_string_pop *l ) { fprintf( f, "\n\tlong_string_pop:\n" ); fprintf( f, "\tGmu_c2 %e\n", l->Gmu_c2 ); fprintf( f, "\tmu %e\n", l->mu ); fprintf( f, "\talpha %e\n", l->alpha ); fprintf( f, "\trho %e\n", l->rho ); fprintf( f, "\td_rho_dt %e\n", l->d_rho_dt); fprintf( f, "\tA %e\n", l->A ); fprintf( f, "\tA_MD %e\n", l->A_MD ); fprintf( f, "\tA_RD %e\n", l->A_RD ); fprintf( f, "\tmean_u2 %e\n", l->mean_u2 ); fprintf( f, "\tmean_u2_RD %e\n", l->mean_u2_RD ); fprintf( f, "\tmean_u2_MD %e\n", l->mean_u2_MD ); } double f(double c1, double v, double Ha) { double fv; fv = c1*v/Ha; return fv; } double g(double c2, double v, double Ha, double ell) { double gv; gv = (1 - SQUARE(v))*(-2*v + c2/(Ha*ell)); return gv; } double hfunc(double c2, double v, double Ha, double ell, double sL) { double hv; hv = c2*sL*v/(2*ell*Ha) + P*ell*v/(2*Ha*sL); return hv; } void rk4( double * y, double * dy, double * yout) { double h; double vtemp, vptemp, elltemp, sLtemp; double k1,k2,k3,k4,l1,l2,l3,l4; double m1,m2,m3,m4,n1,n2,n3,n4; double c1, c2, v, vp, Ha, ell, sL; double newv, newell, newsL, temp; double H, a; c1 = y[0]; c2 = y[1]; H = y[2]; v = y[3]; vp = v*SPEED_OF_LIGHT; ell = y[4]; sL = y[5]; a = y[6]; Ha = H*a; h = LN_A_STEP; k1 = h*f(c1,v,Ha); l1 = h*g(c2,v,Ha,ell); m1 = h*hfunc(c2,v,Ha,ell,sL); elltemp = ell+ k1*0.5; vtemp = v + l1*0.5; sLtemp = sL + m1*0.5; k2 = h*f(c1,vtemp,Ha); l2 = h*g(c2,vtemp,Ha,elltemp); m2 = h*hfunc(c2,vtemp,Ha,elltemp,sLtemp); elltemp = ell + k2*0.5; vtemp = v + l2*0.5; sLtemp = sL + m2*0.5; k3 = h*f(c1,vtemp,Ha); l3 = h*g(c2,vtemp,Ha,elltemp); m3 = h*hfunc(c2,vtemp,Ha,elltemp,sLtemp); elltemp = ell + k3; vtemp = v + l3; sLtemp = sL + m3; k4 = h*f(c1,vtemp,Ha); l4 = h*g(c2,vtemp,Ha,elltemp); m4 = h*hfunc(c2,vtemp,Ha,elltemp,sLtemp); newell = ell + k1/6.0 + k2 / 3.0 + k3 / 3.0 + k4 / 6.0; newv = v + l1/6.0 + l2 / 3.0 + l3 / 3.0 + l4 / 6.0; newsL = sL + m1/6.0 + m2 / 3.0 + m3 / 3.0 + m4 / 6.0; dy[0] = 0; dy[1] = 0; dy[2] = 0; dy[3] = newv - y[1]; dy[4] = newell - y[2]; dy[5] = newsL - y[3]; yout[3] = newv; yout[4] = newell; yout[5] = newsL; }