init_harm_data_old.c

#include "decs.h"

/*

get HARM simulation data from fname

checks for consistency of coordinates in data file with
values of coordinate parameters 

Uses standard HARM data file format

CFG 1 Sept 07

*/

extern double ***bcon;
extern double ***bcov;
extern double ***ucon;
extern double ***ucov;
extern double ***p;
extern double **ne;
extern double **thetae;
extern double **b;

void init_harm_data(char *fname)
{
	FILE *fp;
	double x[4];
	double rp, hp;
	int i, j, k ;

	/* header variables not used except locally */
	double t,tf,cour,DTd,DTl,DTi,dt ;
	int nstep,DTr,dump_cnt,image_cnt,rdump_cnt,lim,failed ;
	double r,h,divb,vmin,vmax,gdet ;
	double J ;
	double dMact, Ladv ;

	fp = fopen(fname, "r");

	if (fp == NULL) {
		fprintf(stderr, "can't open sim data file\n");
		exit(1);
	} else {
		fprintf(stderr,
			"successfully opened %s\n",fname);
	}

	/* get standard HARM header */
	fscanf(fp, "%lf ", &t) ;
	fscanf(fp, "%d ",  &N1) ;
	fscanf(fp, "%d ",  &N2) ;
	fscanf(fp, "%lf ", &startx[1]) ;
	fscanf(fp, "%lf ", &startx[2]) ;
	fscanf(fp, "%lf ", &dx[1]) ;
	fscanf(fp, "%lf ", &dx[2]) ;
	fscanf(fp, "%lf ", &tf) ;
	fscanf(fp, "%d ",  &nstep) ;
	fscanf(fp, "%lf ", &a) ;
	fscanf(fp, "%lf ", &gam) ;
	fscanf(fp, "%lf ", &cour) ;
	fscanf(fp, "%lf ", &DTd) ;
	fscanf(fp, "%lf ", &DTl) ;
	fscanf(fp, "%lf ", &DTi) ;
	fscanf(fp, "%d ",  &DTr) ;
	fscanf(fp, "%d ",  &dump_cnt) ;
	fscanf(fp, "%d ",  &image_cnt) ;
	fscanf(fp, "%d ",  &rdump_cnt) ;
	fscanf(fp, "%lf ", &dt) ;
	fscanf(fp, "%d ",  &lim) ;
	fscanf(fp, "%d ",  &failed) ;
	fscanf(fp, "%lf ", &Rin) ;
	fscanf(fp, "%lf ", &Rout) ;
	fscanf(fp, "%lf ", &hslope) ;
	fscanf(fp, "%lf ", &R0) ;

	//fprintf(stderr,"coordinate parameters: Rin,Rout,hslope,R0,dx[1],dx[2]: %g %g %g %g %g %g\n",
	//	Rin,Rout,hslope,R0,dx[1],dx[2]) ;

	/* nominal non-zero values for axisymmetric simulations */
	startx[0] = 0.;
	startx[3] = 0.;

	stopx[0] = 1. ;
	stopx[1] = startx[1] + N1*dx[1] ;
	stopx[2] = startx[2] + N2*dx[2] ;
	stopx[3] = 2.*M_PI ;

	dx[0] = 1. ;
	dx[3] = 2.*M_PI;

	/* Allocate storage for all model size dependent variables */
	init_storage();

	dMact = 0. ;
	Ladv = 0. ;
	for(k = 0 ; k < N1*N2 ; k++) {
			j = k%N2 ;
			i = (k - j)/N2 ;
			fscanf(fp, "%lf %lf %lf %lf", &x[1], &x[2], &r, &h);

			/* check that we've got the coordinate parameters right */
			bl_coord(x, &rp, &hp);
			if (fabs(rp - r) > 1.e-5 * rp
			    || fabs(hp - h) > 1.e-5) {
				fprintf(stderr, "grid setup error\n");
				fprintf(stderr, "rp,r,hp,h: %g %g %g %g\n",
					rp, r, hp, h);
				fprintf(stderr,
					"edit R0, hslope, compile, and continue\n");
				exit(1);
			}

			fscanf(fp, "%lf %lf %lf %lf %lf %lf %lf %lf",
			       &p[i][j][KRHO], 
			       &p[i][j][UU],
			       &p[i][j][U1], 
			       &p[i][j][U2], 
			       &p[i][j][U3],
			       &p[i][j][B1], 
			       &p[i][j][B2], 
			       &p[i][j][B3]);

			fscanf(fp, "%lf", &divb);

			fscanf(fp, "%lf %lf %lf %lf",
			       &ucon[i][j][0], &ucon[i][j][1],
			       &ucon[i][j][2], &ucon[i][j][3]);
			fscanf(fp, "%lf %lf %lf %lf", &ucov[i][j][0],
			       &ucov[i][j][1], &ucov[i][j][2],
			       &ucov[i][j][3]);
			fscanf(fp, "%lf %lf %lf %lf", &bcon[i][j][0],
			       &bcon[i][j][1], &bcon[i][j][2],
			       &bcon[i][j][3]);
			fscanf(fp, "%lf %lf %lf %lf", &bcov[i][j][0],
			       &bcov[i][j][1], &bcov[i][j][2],
			       &bcov[i][j][3]);
			fscanf(fp, "%lf ", &vmin);
			fscanf(fp, "%lf ", &vmax);
			fscanf(fp, "%lf ", &vmin);
			fscanf(fp, "%lf ", &vmax);
			fscanf(fp, "%lf ", &gdet);

			/* additional stuff: current */
			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf ", &J) ;

			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf ", &J) ;
			fscanf(fp, "%lf\n", &J) ;

			/* check accretion rate */
			if(i <= 20) dMact += gdet * p[i][j][KRHO] * ucon[i][j][1] ;
			if(i >= 20 && i < 40) Ladv += gdet * p[i][j][UU] * ucon[i][j][1] * ucov[i][j][0] ;

		}

	dMact *= dx[3]*dx[2] ;
	dMact /= 21. ;
	Ladv *= dx[3]*dx[2] ;
	Ladv /= 21. ;
	fprintf(stderr,"dMact: %g, Ladv: %g\n",dMact,Ladv) ;


	/* done! */

}