<commit_before>/* 
** profile.c
**
** Program written in order to calculate profile
**
** written by <NAME>
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <malloc.h>
#include <assert.h>
#include <iof.h>

#define ART_LEVEL_ARRAY_LENGTH 10
#define ART_BANNER_LENGTH 45

typedef struct profile_structure {

    int Ntot;
    int Ngas;
    int Ndark;
    int Nstar;
    int Nenctot;
    int Nencgas;
    int Nencdark;
    int Nencstar;
    double ri;
    double rm;
    double ro;
    double vol;
    double Mtot;
    double Mgas;
    double Mdark;
    double Mstar;
    double Menctot;
    double Mencgas;
    double Mencdark;
    double Mencstar;
    double vtot[3];
    double v2tot[6];
    double vgas[3];
    double v2gas[6];
    double vdark[3];
    double v2dark[6];
    double vstar[3];
    double v2star[6];
    double Ltot[3];
    double Lgas[3];
    double Ldark[3];
    double Lstar[3];
    } PS;

typedef struct halo_data {

    int ID;
    int Nbin;
    double rcentre[3];
    double vcentre[3];
    double rbg, Mbg;
    double rcrit, Mcrit;
    double rvcmax, Mrvcmax, vcmax;
    double rmin, rmax;
    PS *ps;
    } HALO_DATA;

typedef struct general_info {

    int positionprecision;
    int gridtype;
    int projectionvariant;
    int centretype;
    int readhalocataloguefile;
    int rmaxfromhalocatalogue;
    int dataformat;
    int halocatalogueformat;
    int verboselevel;
    int Nbin;
    int Nhalo;
    int Nparticleperblockdark, Nparticleinblockdark, Nblockdark;
    int Nparticleperblockstar, Nparticleinblockstar, Nblockstar;
    int Nparticleread;
    int Nrecordread;
    double rhoencbg, rhoenccrit;
    double rmin, rmax;
    double E, Deltabg, Deltacrit;
    double bc[6];
    double binfactor;
    char HaloCatalogueFileName[256];
    char OutputName[256];
    } GI;

typedef struct coordinates {

    double r[3];
    double v[3];
    } COORDINATES;

typedef struct star_properties {

    double mass;
    double initialmass;
    double tform;
    double metallicitySNII;
    double metallicitySNIa;
    } STAR_PROPERTIES;

// => goes to iof
typedef struct cosmological_parameters {

    double ascale0;
    double OmegaM0;
    double OmegaDM0; 
    double OmegaB0;
    double OmegaL0;
    double OmegaK0;
    double OmegaR0;
    double Hubble0;
    double ascale;
    double OmegaM;
    double OmegaDM; 
    double OmegaB;
    double OmegaL;
    double OmegaK;
    double OmegaR;
    double Hubble;
    } COSMOLOGICAL_PARAMETERS;

// => goes to iof
typedef struct unit_system {

    double GNewton;
    double rhocrit0;
    double rhocrit;
    double LBox;
    } UNIT_SYSTEM;


// => goes to iof
typedef struct art_data {

    int doswap;
    int massfromdata;
    int gascontained, darkcontained, starcontained;
    int Lmaxdark;
    int Nparticleperrecord;
    int Nrecord;
    int Ngas, Ndark, Nstar;
    int Ndarklevel[ART_LEVEL_ARRAY_LENGTH];
    int Nstarproperties;
    double shift;
    double toplevelmassdark, toplevelsoftdark, refinementstepdark;
    double massdark[ART_LEVEL_ARRAY_LENGTH];
    double softdark[ART_LEVEL_ARRAY_LENGTH];
    double totalstellarmass, totalstellarinitialmass;
    char HeaderFileName[256], CoordinatesDataFileName[256], StarPropertiesFileName[256];
    char Banner[ART_BANNER_LENGTH];
    ART_HEADER ah;
    FILE *HeaderFile, *CoordinatesDataFile, *StarPropertiesFile;
    } ART_DATA;

typedef struct profile_dark_particle {

    double r[3];
    double v[3];
    double mass;
    } PROFILE_DARK_PARTICLE;

typedef struct profile_star_particle {

    double r[3];
    double v[3];
    double mass;
    double tform;
    double metallicitySNII;
    double metallicitySNIa;
    } PROFILE_STAR_PARTICLE;

void calculate_unit_vectors(double pos[3], double erad[3], double ephi[3], double etheta[3]) {

    double dist;
    double cosphi, sinphi, costheta, sintheta;

    /*
    ** Calculate cosphi & sinphi
    */
    dist = sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
    cosphi = pos[0]/dist;
    sinphi = pos[1]/dist;
    if ((pos[0] == 0) && (pos[1] == 0)) {
        cosphi = 1;
        sinphi = 0;
        }
    /*
    ** Calculate costheta & sintheta
    */
    dist = sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
    costheta = pos[2]/dist;
    sintheta = sqrt(1-costheta*costheta);
    /*
    ** Calculate unit vectors
    */
    erad[0] = sintheta*cosphi;
    erad[1] = sintheta*sinphi;
    erad[2] = costheta;
    ephi[0] = -sinphi;
    ephi[1] = cosphi;
    ephi[2] = 0;
    etheta[0] = -costheta*cosphi;
    etheta[1] = -costheta*sinphi;
    etheta[2] = sintheta;
    }


// => goes to iof
double correct_position(double c, double r, double l) {

    double d;

    d = r-c;
    if (d > 0.5*l) return r - l;
    else if (d < -0.5*l) return r + l;
    else return r;
    }

// => goes to iof
double put_in_box(double r, double lmin, double lmax) {

    double l;

    l = lmax - lmin;
    if (r < lmin) return r + l;
    else if (r > lmax) return r - l;
    else return r;
    }

double Ecosmo(double a, double OmegaM0, double OmegaL0, double OmegaK0) {

    return sqrt(OmegaM0*pow(a,-3.0) + OmegaL0 + OmegaK0*pow(a,-2.0));
    }

void read_art_header(ART_DATA *ad) {

    int i, L;
    int header, trailer;

    ad->doswap = 0;
    if (ad->HeaderFile == NULL) {
	ad->HeaderFile = fopen(ad->HeaderFileName,"r");
	assert(ad->HeaderFile != NULL);
	}
    assert(fread(&header,sizeof(int),1,ad->HeaderFile) == 1);
    if (header != sizeof(ad->Banner)+sizeof(ART_HEADER)) {
	ad->doswap = 1;
	reorder(&header,sizeof(int),1);
	}
    assert(header == sizeof(ad->Banner)+sizeof(ART_HEADER));
    assert(fread(ad->Banner,sizeof(char),sizeof(ad->Banner),ad->HeaderFile) == sizeof(ad->Banner));
    assert(fread(&ad->ah,sizeof(ART_HEADER),1,ad->HeaderFile) == 1);
    if (ad->doswap) reorder(&ad->ah,4,sizeof(ART_HEADER)/4);
    assert(fread(&trailer,sizeof(int),1,ad->HeaderFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    fclose(ad->HeaderFile);
    /*
    ** Set some derived quantities
    */
    ad->Ngas = 0;
    if (ad->starcontained == 1) {
	ad->Lmaxdark = ad->ah.Nspecies-2;
	ad->Ndark = ad->ah.num[ad->Lmaxdark];
	ad->Nstar = ad->ah.num[ad->ah.Nspecies-1] - ad->Ndark;
	}
    else {
	ad->Lmaxdark = ad->ah.Nspecies-1;
	ad->Ndark = ad->ah.num[ad->Lmaxdark];
	ad->Nstar = 0;
	}
    ad->Nparticleperrecord = ad->ah.Nrow*ad->ah.Nrow;
    ad->Nrecord = (ad->Ndark+ad->Nstar+ad->Nparticleperrecord-1)/ad->Nparticleperrecord;
    if (ad->toplevelmassdark == -1) {
	ad->massfromdata = 1;
	ad->toplevelmassdark = ad->ah.mass[ad->Lmaxdark];
	}
    assert(ad->toplevelsoftdark >= 0);
    assert(ad->toplevelmassdark >= 0);
    for (i = 0; i <= ad->Lmaxdark; i++) {
	L = ad->Lmaxdark-i;
	if (i == 0) {
	    ad->Ndarklevel[L] = ad->ah.num[i];
	    }
	else {
	    ad->Ndarklevel[L] = ad->ah.num[i]-ad->ah.num[i-1];
	    }
	if (ad->massfromdata == 1) {
	    ad->massdark[L] = ad->ah.mass[i];
	    }
	}
    }

void read_art_header_star(ART_DATA *ad) {

    int header, trailer;
    int idummy;
    double ddummy;

    if (ad->StarPropertiesFile == NULL) {
	ad->StarPropertiesFile = fopen(ad->StarPropertiesFileName,"r");
	assert(ad->StarPropertiesFile != NULL);
	}

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&idummy,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&idummy,sizeof(int),1);
    assert(idummy == ad->Nstar);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    ad->totalstellarmass = ddummy;
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    ad->totalstellarinitialmass = ddummy;
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    }

void read_art_header_gas(ART_DATA *ad) {

    int header, trailer;
    int idummy;
    double ddummy;

    if (ad->StarPropertiesFile == NULL) {
	ad->StarPropertiesFile = fopen(ad->StarPropertiesFileName,"r");
	assert(ad->StarPropertiesFile != NULL);
	}

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&idummy,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&idummy,sizeof(int),1);
    assert(idummy == ad->Nstar);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);

    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    ad->totalstellarmass = ddummy;
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ddummy,sizeof(double),1);
    ad->totalstellarinitialmass = ddummy;
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    }

void usage(void);
void read_halocatalogue_generic(GI (*), HALO_DATA (**));
void read_halocatalogue_6DFOF(GI (*), HALO_DATA (**));
void initialise_halo_profile (GI (*), HALO_DATA (*));
void read_art_coordinates_record(ART_DATA, COORDINATES (*));
void read_art_star_properties(ART_DATA, int, STAR_PROPERTIES (*));
void put_pdp_in_bins(HALO_DATA (*), PROFILE_DARK_PARTICLE (*), UNIT_SYSTEM, GI);
void put_psp_in_bins(HALO_DATA (*), PROFILE_STAR_PARTICLE (*), UNIT_SYSTEM, GI);
void calculate_halo_properties(HALO_DATA (*), COSMOLOGICAL_PARAMETERS, UNIT_SYSTEM, GI);
void write_output(HALO_DATA (*), GI);

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

    int i, j, k;
    int L = -1;
    int Nparticleread, Nrecordread;
    int Icurrentblockdark, Icurrentblockstar;
    GI gi;
    COSMOLOGICAL_PARAMETERS cp;
    UNIT_SYSTEM us;
    TIPSY_HEADER th;
    GAS_PARTICLE gp;
    DARK_PARTICLE dp;
    STAR_PARTICLE sp;
    GAS_PARTICLE_DPP gpdpp;
    DARK_PARTICLE_DPP dpdpp;
    STAR_PARTICLE_DPP spdpp;
    ART_DATA ad;
    STAR_PROPERTIES starprop;
    COORDINATES *ac = NULL;
    HALO_DATA *hd = NULL;
    PROFILE_DARK_PARTICLE *pdp = NULL;
    PROFILE_STAR_PARTICLE *psp = NULL;
    XDR xdrs;

    /*
    ** Set some default values
    */

    cp.ascale0 = 1;
    cp.OmegaM0 = 0.3;
    cp.OmegaDM0 = 0.26;
    cp.OmegaB0 = 0.04;
    cp.OmegaL0 = 0.7;
    cp.OmegaK0 = 0;
    cp.OmegaR0 = 0;
    cp.Hubble0 = sqrt(8.0*M_PI/3.0);

    us.GNewton = 1;
    us.rhocrit0 = 1;
    us.LBox = 1;

    gi.Deltabg = 200;
    gi.Deltacrit = 0;
    gi.rmin = 0;
    gi.rmax = 1;
    gi.Nbin = 0;
    gi.gridtype = 1;
    gi.binfactor = 5;
    gi.Nhalo = 0;
    gi.Nparticleperblockdark = 10000000;
    gi.Nparticleinblockdark = 0;
    gi.Nblockdark = 0;
    gi.Nparticleperblockstar = 10000000;
    gi.Nparticleinblockstar = 0;
    gi.Nblockstar = 0;

    ad.doswap = 0;
    ad.massfromdata = 0;
    ad.gascontained = 0;
    ad.darkcontained = 0;
    ad.starcontained = 0;
    ad.Ngas = 0;
    ad.Ndark = 0;
    ad.Nstar = 0;
    ad.refinementstepdark = 2;
    ad.toplevelmassdark = -1;
    ad.toplevelsoftdark = 0;
    ad.shift = 0;
    for (i = 0; i < ART_LEVEL_ARRAY_LENGTH; i++) {
	ad.Ndarklevel[i] = 0;
	ad.massdark[i] = 0;
	ad.softdark[i] = 0;
	}
    ad.Nstarproperties = 5;
    ad.HeaderFile = NULL;
    ad.CoordinatesDataFile = NULL;
    ad.StarPropertiesFile = NULL;

    /*
    ** Read in arguments
    */

    i = 1;
    while (i < argc) {
	if (strcmp(argv[i],"-spp") == 0) {
            gi.positionprecision = 0;
            i++;
            }
        else if (strcmp(argv[i],"-dpp") == 0) {
            gi.positionprecision = 1;
            i++;
            }
        else if (strcmp(argv[i],"-lin") == 0) {
	    gi.gridtype = 0;
            i++;
            }
        else if (strcmp(argv[i],"-log") == 0) {
	    gi.gridtype = 1;
            i++;
            }
        else if (strcmp(argv[i],"-v") == 0) {
	    gi.verboselevel = 1;
            i++;
            }
        else if (strcmp(argv[i],"-pv") == 0) {
            i++;
            if (i >= argc) usage();
            gi.projectionvariant = atoi(argv[i]);
            i++;
            }
	else if (strcmp(argv[i],"-rmin") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.rmin = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-rmax") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.rmax = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-binfactor") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.binfactor = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-Nbin") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.Nbin = atof(argv[i]);
	    i++;
	    }
        else if (strcmp(argv[i],"-Delta_bg") == 0) {
            i++;
            if (i >= argc) usage();
            gi.Deltabg = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-Delta_crit") == 0) {
            i++;
            if (i >= argc) usage();
            gi.Deltacrit = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaM0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaM0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaK0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaK0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaL0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaL0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-LBox") == 0) {
            i++;
            if (i >= argc) usage();
            us.LBox = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-rhocrit0") == 0) {
            i++;
            if (i >= argc) usage();
            us.rhocrit0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-G") == 0) {
            i++;
            if (i >= argc) usage();
            us.GNewton = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-com") == 0) {
            gi.centretype = 0;
            i++;
            }
        else if (strcmp(argv[i],"-potmin") == 0) {
            gi.centretype = 1;
            i++;
            }
        else if (strcmp(argv[i],"-denmax") == 0) {
            gi.centretype = 2;
            i++;
            }
        else if (strcmp(argv[i],"-rmaxfromhalocatalogue") == 0) {
            gi.rmaxfromhalocatalogue = 0;
            i++;
            }
        else if (strcmp(argv[i],"-dataformat") == 0) {
            i++;
            if (i >= argc) usage();
	    gi.dataformat = atoi(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-halocatalogueformat") == 0) {
            i++;
            if (i >= argc) usage();
	    gi.halocatalogueformat = atoi(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-halocatalogue") == 0) {
	    gi.readhalocataloguefile = 1;
            i++;
            if (i >= argc) usage();
            strcpy(gi.HaloCatalogueFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-output") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(gi.OutputName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-headerfile") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(ad.HeaderFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-coordinatesdatafile") == 0) {
	    ad.darkcontained = 1;
            i++;
            if (i >= argc) usage();
            strcpy(ad.CoordinatesDataFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-starpropertiesfile") == 0) {
	    ad.starcontained = 1;
            i++;
            if (i >= argc) usage();
            strcpy(ad.StarPropertiesFileName,argv[i]);
            i++;
            }


/*
        else if (strcmp(argv[i],"-gasfile") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(ad.GasFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-darkfile") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(ad.DarkFileName,argv[i]);
            i++;
            }
*/
	else if ((strcmp(argv[i],"-h") == 0) || (strcmp(argv[i],"-help") == 0)) {
	    usage();
	    }
	else {
	    usage();
	    }
	}

    /*
    ** Read header files
    */

    if (gi.dataformat == 0) {
	xdrstdio_create(&xdrs,stdin,XDR_DECODE);
	read_tipsy_standard_header(&xdrs,&th);
	cp.ascale = th.time;
	gi.bc[0] = -0.5*us.LBox;
	gi.bc[1] = -0.5*us.LBox;
	gi.bc[2] = -0.5*us.LBox;
	gi.bc[3] = 0.5*us.LBox;
	gi.bc[4] = 0.5*us.LBox;
	gi.bc[5] = 0.5*us.LBox;
	}
    else if (gi.dataformat == 1) {
	read_art_header(&ad);
	cp.ascale = ad.ah.aunin;
	cp.OmegaM0 = ad.ah.OmM0;
	cp.OmegaB0 = ad.ah.OmB0;
	cp.OmegaDM0 = cp.OmegaM0 - cp.OmegaB0;
	cp.OmegaL0 = ad.ah.OmL0;
	cp.OmegaK0 = ad.ah.OmK0;
	gi.bc[0] = 0;
	gi.bc[1] = 0;
	gi.bc[2] = 0;
	gi.bc[3] = us.LBox;
	gi.bc[4] = us.LBox;
	gi.bc[5] = us.LBox;
	}
    else {
	fprintf(stderr,"Not supported format!\n");
	exit(1);
	}

    fprintf(stderr,"a = %g doswap %d Nparticleperrecord %d %d B %s Ndark %d Lmax %d\n",cp.ascale,ad.doswap,ad.Nparticleperrecord,ad.ah.Nrow,ad.Banner,ad.Ndark,ad.Lmaxdark);
    for (i = 0; i < 10; i++) {
	fprintf(stderr,"i %d Ndarklvel %d massdark %g softdark %g num %d\n",i,ad.Ndarklevel[i],ad.massdark[i],ad.softdark[i],ad.ah.num[i]);
	}

//    exit(1);

    /*
    ** Calculate cosmology relevant stuff => write function for this
    ** Densities in comoving coordinates 
    */

    gi.rhoencbg = gi.Deltabg*cp.OmegaM0*us.rhocrit0;
    gi.E = Ecosmo(cp.ascale,cp.OmegaM0,cp.OmegaL0,cp.OmegaK0);
    cp.OmegaM = cp.OmegaM0/(pow(cp.ascale,3)*gi.E*gi.E);
    us.rhocrit = us.rhocrit0*gi.E*gi.E;
    if (gi.Deltacrit == 0) {
	if (cp.OmegaK0 == 0) {
	    gi.Deltacrit = 178*pow(cp.OmegaM,0.45);
	    }
	else if (cp.OmegaL0 == 0) {
	    gi.Deltacrit = 178*pow(cp.OmegaM,0.3);
	    }
	}
    gi.rhoenccrit = gi.Deltacrit*us.rhocrit*pow(cp.ascale,3);

    fprintf(stderr,"a = %g doswap %d \n",cp.ascale,ad.doswap);

    /*
    ** Read halo catalogue
    */

    if (gi.halocatalogueformat == 0) {
	read_halocatalogue_generic(&gi,&hd);
	}
    else if (gi.halocatalogueformat == 1) {
	read_halocatalogue_6DFOF(&gi,&hd);
	}


    fprintf(stderr,"After halocatalogue\n");
    fprintf(stderr,"file: %s Nhalo %d\n",gi.HaloCatalogueFileName,gi.Nhalo);
    for (i = 0; i < gi.Nhalo; i++) {
	fprintf(stderr,"%d %g %g %g %g %g %g %g %g %d\n",hd[i].ID,hd[i].rcentre[0],hd[i].rcentre[1],hd[i].rcentre[2],hd[i].vcentre[0],hd[i].vcentre[1],hd[i].vcentre[2],hd[i].rmin,hd[i].rmax,hd[i].Nbin);
/*	for(j = 0; j < gi.Nbin+1; j++) {
	    fprintf(stderr,"j %d ri %g ro %g\n",j,hd[i].ps[j].ri,hd[i].ps[j].ro);
	    }
*/
	}

    fprintf(stderr,"Ndark %d Nstar %d\n",ad.Ndark,ad.Nstar);
//    exit(1);

    /*
    ** Harvest data
    */
    if (gi.dataformat == 0) {
	assert(ad.Ngas == 0);
	assert(ad.Nstar == 0);
	assert(gi.positionprecision == 0);
	if (gi.positionprecision == 0) {
	    /*
	    ** Dark Matter
	    */
	    pdp = realloc(pdp,gi.Nparticleperblockdark*sizeof(PROFILE_DARK_PARTICLE));
	    assert(pdp != NULL);
	    Nparticleread = 0;
	    Icurrentblockdark = 0;
	    for (i = 0; i < ad.Ndark; i++) {
		read_tipsy_standard_dark(&xdrs,&dp);
		for (k = 0; k < 3; k++) {
		    pdp[Icurrentblockdark].r[k] = dp.pos[k];
		    pdp[Icurrentblockdark].v[k] = dp.vel[k];
		    }
		pdp[Icurrentblockdark].mass = dp.mass;
		Nparticleread++;
		Icurrentblockdark++;
		if ((Icurrentblockdark == gi.Nparticleperblockdark) || (Nparticleread == ad.Ndark)) {
		    /*
		    ** Block is full or we reached end of dark matter particles
		    */
		    gi.Nparticleinblockdark = Icurrentblockdark;
		    put_pdp_in_bins(hd,pdp,us,gi);
		    Icurrentblockdark = 0;
		    }
		}
	    free(pdp);
	    }
	}
    else if (gi.dataformat == 1) {
	/*
	** Dark Matter and Stars
	*/
	ad.CoordinatesDataFile = fopen(ad.CoordinatesDataFileName,"r");
	assert(ad.CoordinatesDataFile != NULL);
	ac = realloc(ac,ad.Nparticleperrecord*sizeof(COORDINATES));
	assert(ac != NULL);
	pdp = realloc(pdp,gi.Nparticleperblockdark*sizeof(PROFILE_DARK_PARTICLE));
	assert(pdp != NULL);
	psp = realloc(psp,gi.Nparticleperblockstar*sizeof(PROFILE_STAR_PARTICLE));
	assert(psp != NULL);
	Nparticleread = 0;
	Nrecordread = 0;
	Icurrentblockdark = 0;
	Icurrentblockstar = 0;
	/*
	** Read Stars Header
	*/
	read_art_header_star(&ad);
	for (i = 0; i < ad.Nrecord; i++) {
	    read_art_coordinates_record(ad,ac);
	    fprintf(stderr,"i %d Nrecord %d\n",i,ad.Nrecord);
	    for (j = 0; j < ad.Nparticleperrecord; j++) {
		if (Nparticleread < ad.Ndark) {
		    /*
		    ** Dark Matter
		    */
		    for (k = 0; k < 3; k++) {
			pdp[Icurrentblockdark].r[k] = ac[j].r[k] + ad.shift;
			pdp[Icurrentblockdark].v[k] = ac[j].v[k];
			}
		    for (k = ad.Lmaxdark; k >=0; k--) {
			if (ad.ah.num[k] >= Nparticleread) L = ad.Lmaxdark-k;
			}
		    pdp[Icurrentblockdark].mass = ad.massdark[L];
		    Nparticleread++;
		    Icurrentblockdark++;
		    if ((Icurrentblockdark == gi.Nparticleperblockdark) || (Nparticleread == ad.Ndark)) {
			/*
			** Block is full or we reached end of dark matter particles
			*/
			gi.Nparticleinblockdark = Icurrentblockdark;<commit_msg>Security update. Implemented sfc library stuff now.<commit_after>/* 
** profile.c
**
** Program written in order to calculate profile
**
** written by <NAME>
*/

#include <stdio.h>
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <malloc.h>
#include <assert.h>
#include <iof.h>
#include <sfc.h>

#define ART_NUMBER_GAS_LEVELS 30
#define ART_NUMBER_DARK_LEVELS 10
#define ART_BANNER_LENGTH 45

typedef struct profile_structure {

    int Ntot;
    int Ngas;
    int Ndark;
    int Nstar;
    int Nenctot;
    int Nencgas;
    int Nencdark;
    int Nencstar;
    double ri;
    double rm;
    double ro;
    double vol;
    double Mtot;
    double Mgas;
    double Mdark;
    double Mstar;
    double Menctot;
    double Mencgas;
    double Mencdark;
    double Mencstar;
    double vtot[3];
    double v2tot[6];
    double vgas[3];
    double v2gas[6];
    double vdark[3];
    double v2dark[6];
    double vstar[3];
    double v2star[6];
    double Ltot[3];
    double Lgas[3];
    double Ldark[3];
    double Lstar[3];
    } PS;

typedef struct halo_data {

    int ID;
    int Nbin;
    double rcentre[3];
    double vcentre[3];
    double rbg, Mbg;
    double rcrit, Mcrit;
    double rvcmax, Mrvcmax, vcmax;
    double rmin, rmax;
    PS *ps;
    } HALO_DATA;

typedef struct general_info {

    int positionprecision;
    int gridtype;
    int projectionvariant;
    int centretype;
    int readhalocataloguefile;
    int rmaxfromhalocatalogue;
    int dataformat;
    int halocatalogueformat;
    int verboselevel;
    int Nbin;
    int Nhalo;
    int Nparticleperblockdark, Nparticleinblockdark, Nblockdark;
    int Nparticleperblockstar, Nparticleinblockstar, Nblockstar;
    int Nparticleread;
    int Nrecordread;
    double rhoencbg, rhoenccrit;
    double rmin, rmax;
    double E, Deltabg, Deltacrit;
    double bc[6];
    double binfactor;
    char HaloCatalogueFileName[256];
    char OutputName[256];
    } GI;

typedef struct coordinates {

    double r[3];
    double v[3];
    } COORDINATES;

typedef struct star_properties {

    double mass;
    double initialmass;
    double tform;
    double metallicitySNII;
    double metallicitySNIa;
    } STAR_PROPERTIES;

// => goes to iof
typedef struct cosmological_parameters {

    double ascale0;
    double OmegaM0;
    double OmegaDM0; 
    double OmegaB0;
    double OmegaL0;
    double OmegaK0;
    double OmegaR0;
    double Hubble0;
    double ascale;
    double OmegaM;
    double OmegaDM; 
    double OmegaB;
    double OmegaL;
    double OmegaK;
    double OmegaR;
    double Hubble;
    } COSMOLOGICAL_PARAMETERS;

// => goes to iof
typedef struct unit_system {

    double GNewton;
    double rhocrit0;
    double rhocrit;
    double LBox;
    } UNIT_SYSTEM;


// => goes to iof
typedef struct art_data {

    /*
    ** from general header
    */
    char Banner[ART_BANNER_LENGTH];
    ART_HEADER ah;
    /*
    ** from star properties file
    */
    double totalstellarmass, totalstellarinitialmass;
    /*
    ** from gas file
    */
    float refinementvolumemin[3], refinementvolumemax[3];
    float starformationvolumemin[3], starformationvolumemax[3];
    SFC_INFO sfci;
    /*
    ** new and derived stuff to get data better organised
    */
    int doswap;
    int massfromdata;
    int gascontained, darkcontained, starcontained;
    int Ndim;
    int Lmaxdark;
    int Lmingas, Lmaxgas, Nlevelgas;
    int Nparticleperrecord;
    int Nrecord;
    int Ngas, Ndark, Nstar;
    int Ndarklevel[ART_NUMBER_DARK_LEVELS];
    int Nstarproperties;
    int Nhydroproperties, Notherproperties;
    int Nrtchemspecies, Nchemspecies, Nrtdiskvars;
    int GRAVITY, HYDRO, STARFORM, ADVECT_SPECIES, ENRICH, ENRICH_SNIa, RADIATIVE_TRANSFER, ELECTRON_ION_NONEQUILIBRIUM;
    long Ncell[ART_NUMBER_GAS_LEVELS];
    double shift;
    double toplevelmassdark, toplevelsoftdark, refinementstepdark;
    double massdark[ART_NUMBER_DARK_LEVELS];
    double softdark[ART_NUMBER_DARK_LEVELS];
    char HeaderFileName[256], CoordinatesDataFileName[256], StarPropertiesFileName[256], GasFileName[256];
    FILE *HeaderFile, *CoordinatesDataFile, *StarPropertiesFile, *GasFile;
    } ART_DATA;

typedef struct profile_dark_particle {

    double r[3];
    double v[3];
    double mass;
    } PROFILE_DARK_PARTICLE;

typedef struct profile_star_particle {

    double r[3];
    double v[3];
    double mass;
    double tform;
    double metallicitySNII;
    double metallicitySNIa;
    } PROFILE_STAR_PARTICLE;

void calculate_unit_vectors(double pos[3], double erad[3], double ephi[3], double etheta[3]) {

    double dist;
    double cosphi, sinphi, costheta, sintheta;

    /*
    ** Calculate cosphi & sinphi
    */
    dist = sqrt(pos[0]*pos[0]+pos[1]*pos[1]);
    cosphi = pos[0]/dist;
    sinphi = pos[1]/dist;
    if ((pos[0] == 0) && (pos[1] == 0)) {
        cosphi = 1;
        sinphi = 0;
        }
    /*
    ** Calculate costheta & sintheta
    */
    dist = sqrt(pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2]);
    costheta = pos[2]/dist;
    sintheta = sqrt(1-costheta*costheta);
    /*
    ** Calculate unit vectors
    */
    erad[0] = sintheta*cosphi;
    erad[1] = sintheta*sinphi;
    erad[2] = costheta;
    ephi[0] = -sinphi;
    ephi[1] = cosphi;
    ephi[2] = 0;
    etheta[0] = -costheta*cosphi;
    etheta[1] = -costheta*sinphi;
    etheta[2] = sintheta;
    }


// => goes to iof
double correct_position(double c, double r, double l) {

    double d;

    d = r-c;
    if (d > 0.5*l) return r - l;
    else if (d < -0.5*l) return r + l;
    else return r;
    }

// => goes to iof
double put_in_box(double r, double lmin, double lmax) {

    double l;

    l = lmax - lmin;
    if (r < lmin) return r + l;
    else if (r > lmax) return r - l;
    else return r;
    }

double Ecosmo(double a, double OmegaM0, double OmegaL0, double OmegaK0) {

    return sqrt(OmegaM0*pow(a,-3.0) + OmegaL0 + OmegaK0*pow(a,-2.0));
    }

void read_art_header(ART_DATA *ad) {

    int i, L;
    int header, trailer;

    ad->doswap = 0;
    if (ad->HeaderFile == NULL) {
	ad->HeaderFile = fopen(ad->HeaderFileName,"r");
	assert(ad->HeaderFile != NULL);
	}
    assert(fread(&header,sizeof(int),1,ad->HeaderFile) == 1);
    if (header != sizeof(ad->Banner)+sizeof(ART_HEADER)) {
	ad->doswap = 1;
	reorder(&header,sizeof(int),1);
	}
    assert(header == sizeof(ad->Banner)+sizeof(ART_HEADER));
    assert(fread(ad->Banner,sizeof(char),sizeof(ad->Banner),ad->HeaderFile) == sizeof(ad->Banner));
    assert(fread(&ad->ah,sizeof(ART_HEADER),1,ad->HeaderFile) == 1);
    if (ad->doswap) reorder(&ad->ah,4,sizeof(ART_HEADER)/4);
    assert(fread(&trailer,sizeof(int),1,ad->HeaderFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    fclose(ad->HeaderFile);
    /*
    ** Set some derived quantities
    */
    ad->Ngas = 0;
    if (ad->starcontained == 1) {
	ad->Lmaxdark = ad->ah.Nspecies-2;
	ad->Ndark = ad->ah.num[ad->Lmaxdark];
	ad->Nstar = ad->ah.num[ad->ah.Nspecies-1] - ad->Ndark;
	}
    else {
	ad->Lmaxdark = ad->ah.Nspecies-1;
	ad->Ndark = ad->ah.num[ad->Lmaxdark];
	ad->Nstar = 0;
	}
    ad->Nparticleperrecord = ad->ah.Nrow*ad->ah.Nrow;
    ad->Nrecord = (ad->Ndark+ad->Nstar+ad->Nparticleperrecord-1)/ad->Nparticleperrecord;
    if (ad->toplevelmassdark == -1) {
	ad->massfromdata = 1;
	ad->toplevelmassdark = ad->ah.mass[ad->Lmaxdark];
	}
    assert(ad->toplevelsoftdark >= 0);
    assert(ad->toplevelmassdark >= 0);
    for (i = 0; i <= ad->Lmaxdark; i++) {
	L = ad->Lmaxdark-i;
	if (i == 0) {
	    ad->Ndarklevel[L] = ad->ah.num[i];
	    }
	else {
	    ad->Ndarklevel[L] = ad->ah.num[i]-ad->ah.num[i-1];
	    }
	if (ad->massfromdata == 1) {
	    ad->massdark[L] = ad->ah.mass[i];
	    }
	}
    }

void read_art_header_star(ART_DATA *ad) {

    int header, trailer;
    int idummy;
    double ddummy;

    if (ad->StarPropertiesFile == NULL) {
	ad->StarPropertiesFile = fopen(ad->StarPropertiesFileName,"r");
	assert(ad->StarPropertiesFile != NULL);
	}
    /*
    ** st, sa
    */
    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    assert(fread(&ddummy,sizeof(double),1,ad->StarPropertiesFile) == 1);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** Nstar
    */
    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&idummy,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&idummy,sizeof(int),1);
    assert(idummy == ad->Nstar);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** total stellar mass, total stellar initial mass
    */
    assert(fread(&header,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ad->totalstellarmass,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ad->totalstellarmass,sizeof(double),1);
    assert(fread(&ad->totalstellarinitialmass,sizeof(double),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&ad->totalstellarinitialmass,sizeof(double),1);
    assert(fread(&trailer,sizeof(int),1,ad->StarPropertiesFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    }

void read_art_header_gas(ART_DATA *ad) {

    int header, trailer, i;
    int idummy;
    float fdummy;
    double ddummy;
    char cdummy[256];

    if (ad->GasFile == NULL) {
	ad->GasFile = fopen(ad->GasFileName,"r");
	assert(ad->GasFile != NULL);
	}
    /*
    ** Jobname
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&cdummy,sizeof(char),256,ad->GasFile) == 256);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** istep, t, dt, adum, ainit 
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&idummy,sizeof(int),1,ad->GasFile) == 1);
    assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
    assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** boxh, Om0, Oml0, Omb0, hubble
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&fdummy,sizeof(float),1,ad->GasFile) == 1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** Nextra (some old crap - should be zero!)
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&idummy,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&idummy,sizeof(int),1);
    assert(idummy == 0);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** extra
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** lextra
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** Minlevel, Maxlevel
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ad->Lmingas,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&ad->Lmingas,sizeof(int),1);
    assert(fread(&ad->Lmaxgas,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&ad->Lmaxgas,sizeof(int),1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    ad->Nlevelgas = ad->Lmaxgas-ad->Lmingas+1;
    /*
    ** tl
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    for (i = 0; i < ad->Nlevelgas; i++) {
	assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
	}
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** dtl
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    for (i = 0; i < ad->Nlevelgas; i++) {
	assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
	}
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /*
    ** tl_old
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    for (i = 0; i < ad->Nlevelgas; i++) {
	assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
	}
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** dtl_old
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    for (i = 0; i < ad->Nlevelgas; i++) {
	assert(fread(&ddummy,sizeof(double),1,ad->GasFile) == 1);
	}
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** iSO (sweep direction for flux solver)
    */
    if (ad->HYDRO) {
	assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
	if (ad->doswap) reorder(&header,sizeof(int),1);
	for (i = 0; i < ad->Nlevelgas; i++) {
	    assert(fread(&idummy,sizeof(int),1,ad->GasFile) == 1);
	    }
	assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
	if (ad->doswap) reorder(&trailer,sizeof(int),1);
	assert(header == trailer);
	}
    /*
    ** space filling curve (SFC) order
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ad->sfci.sfc_order,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&ad->sfci.sfc_order,sizeof(int),1);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** refinement volume 
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ad->refinementvolumemin,sizeof(float),ad->Ndim,ad->GasFile) == ad->Ndim);
    if (ad->doswap) reorder(&ad->refinementvolumemin,sizeof(float),ad->Ndim);
    assert(fread(&ad->refinementvolumemax,sizeof(float),ad->Ndim,ad->GasFile) == ad->Ndim);
    if (ad->doswap) reorder(&ad->refinementvolumemax,sizeof(float),ad->Ndim);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    /* 
    ** star formation volume 
    */
    if (ad->STARFORM) {
	assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
	if (ad->doswap) reorder(&header,sizeof(int),1);
	assert(fread(&ad->starformationvolumemin,sizeof(float),ad->Ndim,ad->GasFile) == ad->Ndim);
	if (ad->doswap) reorder(&ad->starformationvolumemin,sizeof(float),ad->Ndim);
	assert(fread(&ad->starformationvolumemax,sizeof(float),ad->Ndim,ad->GasFile) == ad->Ndim);
	if (ad->doswap) reorder(&ad->starformationvolumemax,sizeof(float),ad->Ndim);
	assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
	if (ad->doswap) reorder(&trailer,sizeof(int),1);
	assert(header == trailer);
	}
    /*
    ** Ncell[0]
    */
    assert(fread(&header,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&header,sizeof(int),1);
    assert(fread(&ad->Ncell[0],sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&ad->Ncell[0],sizeof(int),1);
    assert(ad->Ncell[0] == ad->ah.Ngrid*ad->ah.Ngrid*ad->ah.Ngrid);
    assert(fread(&trailer,sizeof(int),1,ad->GasFile) == 1);
    if (ad->doswap) reorder(&trailer,sizeof(int),1);
    assert(header == trailer);
    }

void usage(void);
void read_halocatalogue_generic(GI (*), HALO_DATA (**));
void read_halocatalogue_6DFOF(GI (*), HALO_DATA (**));
void initialise_halo_profile (GI (*), HALO_DATA (*));
void read_art_coordinates_record(ART_DATA, COORDINATES (*));
void read_art_star_properties(ART_DATA, int, STAR_PROPERTIES (*));
void put_pdp_in_bins(HALO_DATA (*), PROFILE_DARK_PARTICLE (*), UNIT_SYSTEM, GI);
void put_psp_in_bins(HALO_DATA (*), PROFILE_STAR_PARTICLE (*), UNIT_SYSTEM, GI);
void calculate_halo_properties(HALO_DATA (*), COSMOLOGICAL_PARAMETERS, UNIT_SYSTEM, GI);
void write_output(HALO_DATA (*), GI);

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

    int i, j, k;
    int L = -1;
    int Nparticleread, Nrecordread;
    int Icurrentblockdark, Icurrentblockstar;
    GI gi;
    COSMOLOGICAL_PARAMETERS cp;
    UNIT_SYSTEM us;
    TIPSY_HEADER th;
    GAS_PARTICLE gp;
    DARK_PARTICLE dp;
    STAR_PARTICLE sp;
    GAS_PARTICLE_DPP gpdpp;
    DARK_PARTICLE_DPP dpdpp;
    STAR_PARTICLE_DPP spdpp;
    ART_DATA ad;
    STAR_PROPERTIES starprop;
    COORDINATES *ac = NULL;
    HALO_DATA *hd = NULL;
    PROFILE_DARK_PARTICLE *pdp = NULL;
    PROFILE_STAR_PARTICLE *psp = NULL;
    XDR xdrs;

    /*
    ** Set some default values
    */

    cp.ascale0 = 1;
    cp.OmegaM0 = 0.3;
    cp.OmegaDM0 = 0.26;
    cp.OmegaB0 = 0.04;
    cp.OmegaL0 = 0.7;
    cp.OmegaK0 = 0;
    cp.OmegaR0 = 0;
    cp.Hubble0 = sqrt(8.0*M_PI/3.0);

    us.GNewton = 1;
    us.rhocrit0 = 1;
    us.LBox = 1;

    gi.Deltabg = 200;
    gi.Deltacrit = 0;
    gi.rmin = 0;
    gi.rmax = 1;
    gi.Nbin = 0;
    gi.gridtype = 1;
    gi.binfactor = 5;
    gi.Nhalo = 0;
    gi.Nparticleperblockdark = 10000000;
    gi.Nparticleinblockdark = 0;
    gi.Nblockdark = 0;
    gi.Nparticleperblockstar = 10000000;
    gi.Nparticleinblockstar = 0;
    gi.Nblockstar = 0;

    ad.doswap = 0;
    ad.massfromdata = 0;
    ad.gascontained = 0;
    ad.darkcontained = 0;
    ad.starcontained = 0;
    ad.Ndim = 3;
    ad.Ngas = 0;
    ad.Ndark = 0;
    ad.Nstar = 0;
    ad.refinementstepdark = 2;
    ad.toplevelmassdark = -1;
    ad.toplevelsoftdark = 0;
    ad.shift = 0;
    for (i = 0; i < ART_NUMBER_DARK_LEVELS; i++) {
	ad.Ndarklevel[i] = 0;
	ad.massdark[i] = 0;
	ad.softdark[i] = 0;
	}
    for (i = 0; i < ART_NUMBER_GAS_LEVELS; i++) {
	ad.Ncell[i] = 0;
	}
    ad.Nstarproperties = 5;
    ad.HeaderFile = NULL;
    ad.CoordinatesDataFile = NULL;
    ad.StarPropertiesFile = NULL;
    ad.GRAVITY = 1;
    ad.HYDRO = 1;
    ad.STARFORM = 1;
    ad.ADVECT_SPECIES = 1;
    ad.ENRICH = 1;
    ad.ENRICH_SNIa = 1;
    ad.RADIATIVE_TRANSFER = 1;
    ad.ELECTRON_ION_NONEQUILIBRIUM = 0;

    /*
    ** Read in arguments
    */

    i = 1;
    while (i < argc) {
	if (strcmp(argv[i],"-spp") == 0) {
            gi.positionprecision = 0;
            i++;
            }
        else if (strcmp(argv[i],"-dpp") == 0) {
            gi.positionprecision = 1;
            i++;
            }
        else if (strcmp(argv[i],"-lin") == 0) {
	    gi.gridtype = 0;
            i++;
            }
        else if (strcmp(argv[i],"-log") == 0) {
	    gi.gridtype = 1;
            i++;
            }
        else if (strcmp(argv[i],"-v") == 0) {
	    gi.verboselevel = 1;
            i++;
            }
        else if (strcmp(argv[i],"-pv") == 0) {
            i++;
            if (i >= argc) usage();
            gi.projectionvariant = atoi(argv[i]);
            i++;
            }
	else if (strcmp(argv[i],"-rmin") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.rmin = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-rmax") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.rmax = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-binfactor") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.binfactor = atof(argv[i]);
	    i++;
	    }
	else if (strcmp(argv[i],"-Nbin") == 0) {
	    i++;
            if (i >= argc) usage();
	    gi.Nbin = atof(argv[i]);
	    i++;
	    }
        else if (strcmp(argv[i],"-Delta_bg") == 0) {
            i++;
            if (i >= argc) usage();
            gi.Deltabg = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-Delta_crit") == 0) {
            i++;
            if (i >= argc) usage();
            gi.Deltacrit = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaM0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaM0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaK0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaK0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-OmegaL0") == 0) {
            i++;
            if (i >= argc) usage();
            cp.OmegaL0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-LBox") == 0) {
            i++;
            if (i >= argc) usage();
            us.LBox = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-rhocrit0") == 0) {
            i++;
            if (i >= argc) usage();
            us.rhocrit0 = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-G") == 0) {
            i++;
            if (i >= argc) usage();
            us.GNewton = atof(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-com") == 0) {
            gi.centretype = 0;
            i++;
            }
        else if (strcmp(argv[i],"-potmin") == 0) {
            gi.centretype = 1;
            i++;
            }
        else if (strcmp(argv[i],"-denmax") == 0) {
            gi.centretype = 2;
            i++;
            }
        else if (strcmp(argv[i],"-rmaxfromhalocatalogue") == 0) {
            gi.rmaxfromhalocatalogue = 0;
            i++;
            }
        else if (strcmp(argv[i],"-dataformat") == 0) {
            i++;
            if (i >= argc) usage();
	    gi.dataformat = atoi(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-halocatalogueformat") == 0) {
            i++;
            if (i >= argc) usage();
	    gi.halocatalogueformat = atoi(argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-halocatalogue") == 0) {
	    gi.readhalocataloguefile = 1;
            i++;
            if (i >= argc) usage();
            strcpy(gi.HaloCatalogueFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-output") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(gi.OutputName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-headerfile") == 0) {
            i++;
            if (i >= argc) usage();
            strcpy(ad.HeaderFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-coordinatesdatafile") == 0) {
	    ad.darkcontained = 1;
            i++;
            if (i >= argc) usage();
            strcpy(ad.CoordinatesDataFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-starpropertiesfile") == 0) {
	    ad.starcontained = 1;
            i++;
            if (i >= argc) usage();
            strcpy(ad.StarPropertiesFileName,argv[i]);
            i++;
            }
        else if (strcmp(argv[i],"-gasfile") == 0) {
	    ad.gascontained = 1;
            i++;
            if (i >= argc) usage();
            strcpy(ad.GasFileName,argv[i]);
            i++;
            }
	else if ((strcmp(argv[i],"-h") == 0) || (strcmp(argv[i],"-help") == 0)) {
	    usage();
	    }
	else {
	    usage();
	    }
	}

    /*
    ** Read header files
    */

    if (gi.dataformat == 0) {
	xdrstdio_create(&xdrs,stdin,XDR_DECODE);
	read_tipsy_standard_header(&xdrs,&th);
	cp.ascale = th.time;
	gi.bc[0] = -0.5*us.LBox;
	gi.bc[1] = -0.5*us.LBox;
	gi.bc[2] = -0.5*us.LBox;
	gi.bc[3] = 0.5*us.LBox;
	gi.bc[4] = 0.5*us.LBox;
	gi.bc[5] = 0.5*us.LBox;
	}
    else if (gi.dataformat == 1) {
	read_art_header(&ad);
	cp.ascale = ad.ah.aunin;
	cp.OmegaM0 = ad.ah.OmM0;
	cp.OmegaB0 = ad.ah.OmB0;
	cp.OmegaDM0 = cp.OmegaM0 - cp.OmegaB0;
	cp.OmegaL0 = ad.ah.OmL0;
	cp.OmegaK0 = ad.ah.OmK0;
	gi.bc[0] = 0;
	gi.bc[1] = 0;
	gi.bc[2] = 0;
	gi.bc[3] = us.LBox;
	gi.bc[4] = us.LBox;
	gi.bc[5] = us.LBox;
	}
    else {
	fprintf(stderr,"Not supported format!\n");
	exit(1);
	}

    fprintf(stderr,"a = %g doswap %d Nparticleperrecord %d %d B %s Ndark %d Lmax %d\n",cp.ascale,ad.doswap,ad.Nparticleperrecord,ad.ah.Nrow,ad.Banner,ad.Ndark,ad.Lmaxdark);
    for (i = 0; i < 10; i++) {
	fprintf(stderr,"i %d Ndarklvel %d massdark %g softdark %g num %d\n",i,ad.Ndarklevel[i],ad.massdark[i],ad.softdark[i],ad.ah.num[i]);
	}

//    exit(1);

    /*
    ** Calculate cosmology relevant stuff => write function for this
    ** Densities in comoving coordinates 
    */

    gi.rhoencbg = gi.Deltabg*cp.OmegaM0*us.rhocrit0;
    gi.E = Ecosmo(cp.ascale,cp.OmegaM0,cp.OmegaL0,cp.OmegaK0);
    cp.OmegaM = cp.OmegaM0/(pow(cp.ascale,3)*gi.E*gi.E);
    us.rhocrit = us.rhocrit0*gi.E*gi.E;
    if (gi.Deltacrit == 0) {
	if (cp.OmegaK0 == 0) {
	    gi.Deltacrit = 178*pow(cp.OmegaM,0.45);
	    }
	else if (cp.OmegaL0 == 0) {
	    gi.Deltacrit = 178*pow(cp.OmegaM,0.3);
	    }
	}
    gi.rhoenccrit = gi.Deltacrit*us.rhocrit*pow(cp.ascale,3);

    fprintf(stderr,"a = %g doswap %d \n",cp.ascale,ad.doswap);

    /*
    ** Read halo catalogue
    */

    if (gi.halocatalogueformat == 0) {
	read_halocatalogue_generic(&gi,&hd);
	}
    else if (gi.halocatalogueformat == 1) {
	read_halocatalogue_6DFOF(&gi,&hd);
	}


    fprintf(stderr,"After halocatalogue\n");
    fprintf(stderr,"file: %s Nhalo %d\n",gi.HaloCatalogueFileName,gi.Nhalo);
    for (i = 0; i < gi.Nhalo; i++) {
	fprintf(stderr,"%d %g %g %g %g %g %g %g %g %d\n",hd[i].ID,hd[i].rcentre[0],hd[i].rcentre[1],hd[i].rcentre[2],hd[i].vcentre[0],hd[i].vcentre[1],hd[i].vcentre[2],hd[i].rmin,hd[i].rmax,hd[i].Nbin);
/*	for(j = 0; j < gi.Nbin+1; j++) {
	    fprintf(stderr,"j %d ri %g ro %g\n",j,hd[i].ps[j].ri,hd[i].ps[j].ro);
	    }
*/
	}

    fprintf(stderr,"Ndark %d Nstar %d\n",ad.Ndark,ad.Nstar);
//    exit(1);

    /*
    ** Harvest data
    */
    if (gi.dataformat == 0) {
	assert(ad.Ngas == 0);
	assert(ad.Nstar == 0);
	assert(gi.positionprecision == 0);
	if (gi.positionprecision == 0) {
	    /*
	    ** Dark Matter
	    */
	    pdp = realloc(pdp,gi.Nparticleperblockdark*sizeof(PROFILE_DARK_PARTICLE));
	    assert(pdp != NULL);
	    Nparticleread = 0;
	    Icurrentblockdark = 0;
	    for (i = 0; i < ad.Ndark; i++) {
		read_tipsy_standard_dark(&xdrs,&dp);
		for (k = 0; k < 3; k++) {
		    pdp[Icurrentblockdark].r[k] = dp.pos[k];
		    pdp[Icurrentblockdark].v[k] = dp.vel[k];
		    }
		pdp[Icurrentblockdark].mass = dp.mass;
		Nparticleread++;
		Icurrentblockdark++;
		if ((Icurrentblockdark == gi.Nparticleperblockdark) || (Nparticleread == ad.Ndark)) {
		    /*
		    ** Block is full or we reached end of dark matter particles
		    */
		    gi.Nparticleinblockdark = Icurrentblockdark;
		    put_pdp_in_bins(hd,pdp,us,gi);
		    Icurrentblockdark = 0;
		    }
		}
	    free(pdp);
	    }
	}
    else if (gi.dataformat == 1) {
	/*
	** Gas
	*/
	ad.GasFile = fopen(ad.GasFileName,"r");
	assert(ad.GasFile != NULL);
	read_art_header_gas(&ad);

	ad.Nhydroproperties = 0;
	ad.Notherproperties = 0;
	ad.Nrtchemspecies = 0;
	ad.Nchemspecies = 0;
	ad.Nrtdiskvars = 0;
	if (ad.HYDRO) {
	    if (ad.ADVECT_SPECIES) {
		if (ad.RADIATIVE_TRANSFER) ad.Nrtchemspecies = 6;
		else ad.Nrtchemspecies = 0;
		if (ad.ENRICH) {
		    if (ad.ENRICH_SNIa) ad.Nchemspecies = ad.Nrtchemspecies + 2;
		    else ad.Nchemspecies = ad.Nrtchemspecies + 1;
		    }
		else ad.Nchemspecies = ad.Nrtchemspecies;
		}
	    else ad.Nchemspecies = 0;
	    if (ad.ELECTRON_ION_NONEQUILIBRIUM) ad.Nhydroproperties = 6 + ad.Ndim + ad.Nchemspecies;
	    else ad.Nhydroproperties = 5 + ad.Ndim + ad.Nchemspecies;
	    }
	if (ad.RADIATIVE_TRANSFER) ad.Nrtdiskvars = 6;
	else ad.Nrtdiskvars = 0;
	if (ad.GRAVITY) ad.Notherproperties++;
	if (ad.HYDRO) ad.Notherproperties++;
	if (ad.RADIATIVE_TRANSFER) ad.Notherproperties += ad.Nrtdiskvars;

	fprintf(stderr,"B: Ndim %d Ngrid %d nBits %d nBitsPerDim %d\n",ad.Ndim,ad.ah.Ngrid,ad.sfci.nBits,ad.sfci.nBitsPerDim);

	init_sfc(&ad.sfci,ad.Ndim,ad.ah.Ngrid);

	fprintf(stderr,"A: Ndim %d Ngrid %d nBits %d nBitsPerDim %d\n",ad.Ndim,ad.ah.Ngrid,ad.sfci.nBits,ad.sfci.nBitsPerDim);

	exit(1);

	int header,trailer;
	int *cellrefinedbuffer = NULL;
	float cellhydroproperties[ad.Nhydroproperties];
	float cellotherproperties[ad.Notherproperties];

	fprintf(stderr,"Nhydroproperties %d Notherproperties %d sfc_order %d\n",
		ad.Nhydroproperties,ad.Notherproperties,ad.sfci.sfc_order);

	assert(fread(&header,sizeof(int),1,ad.GasFile) == 1);
	if (ad.doswap) reorder(&header,sizeof(int),1);
	if (ad.Lmingas == 0) assert(ad.Ncell[0] == header/sizeof(int));
	ad.Ncell[ad.Lmingas] = header/(sizeof(int));
	cellrefinedbuffer = realloc(cellrefinedbuffer,ad.Ncell[ad.Lmingas]*sizeof(int));
	assert(fread(cellrefinedbuffer,sizeof(int),ad.Ncell[ad.Lmingas],ad.GasFile) == ad.Ncell[ad.Lmingas]);
	if (ad.doswap) reorder(&cellrefinedbuffer,sizeof(int),ad.Ncell[ad.Lmingas]);
	assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	if (ad.doswap) reorder(&trailer,sizeof(int),1);
	assert(header == trailer);

	fprintf(stderr,"Ncell[0] %lu =? %lu h %d t %d Lmingas %d Lgasmax %d\n",
		ad.Ncell[0],header/sizeof(int),header,trailer,ad.Lmingas,ad.Lmaxgas);
/*
	for (i = 0; i < ad.Ncell[0]; i++) {
	    if (cellrefinedbuffer[i] > 1) {
		fprintf(stderr,"i %d crb %d\n",i,cellrefinedbuffer[i]);
		}
	    }
*/


	/*
	** min level hydro properties
	*/
	fread(&header,sizeof(int),1,ad.GasFile);
	if (ad.doswap) reorder(&header,sizeof(int),1);
	assert(ad.Nhydroproperties == header/(sizeof(float)*ad.Ncell[ad.Lmingas]));
	for (i = 0; i < ad.Ncell[ad.Lmingas]; i++) {
	    assert(fread(cellhydroproperties,sizeof(float),ad.Nhydroproperties,ad.GasFile) == ad.Nhydroproperties);
	    if (ad.doswap) reorder(cellhydroproperties,sizeof(float),ad.Nhydroproperties);
	    }
	assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	if (ad.doswap) reorder(&trailer,sizeof(int),1);
	assert(header == trailer);
//	fprintf(stderr,"h %d t %d\n",header,trailer);
	/*
	** min level other properties
	*/
	if (ad.GRAVITY || ad.RADIATIVE_TRANSFER) {
	    assert(fread(&header,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&header,sizeof(int),1);
	    assert(ad.Notherproperties == header/(sizeof(float)*ad.Ncell[ad.Lmingas]));
	    for (i = 0; i < ad.Ncell[ad.Lmingas]; i++) {
		assert(fread(cellotherproperties,sizeof(float),ad.Notherproperties,ad.GasFile) == ad.Notherproperties);
		if (ad.doswap) reorder(cellotherproperties,sizeof(float),ad.Notherproperties);
		}
	    assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&trailer,sizeof(int),1);
	    assert(header == trailer);
//	    fprintf(stderr,"h %d t %d\n",header,trailer);
	    }

	for (i = ad.Lmingas+1; i <= ad.Lmaxgas; i++) {
	    /*
	    ** number of cells
	    */
	    assert(fread(&header,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&header,sizeof(int),1);
	    assert(fread(&ad.Ncell[i],sizeof(long),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&ad.Ncell[i],sizeof(long),1);
	    assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&trailer,sizeof(int),1);
//	    fprintf(stderr,"h %d t %d\n",header,trailer);
	    assert(header == trailer);
	    /*
	    ** cellrefinedbuffer
	    */
	    cellrefinedbuffer = realloc(cellrefinedbuffer,ad.Ncell[i]*sizeof(int));
	    assert(fread(&header,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&header,sizeof(int),1);
	    assert(fread(cellrefinedbuffer,sizeof(int),ad.Ncell[i],ad.GasFile) == ad.Ncell[i]);
	    if (ad.doswap) reorder(&cellrefinedbuffer,sizeof(int),ad.Ncell[i]);
	    assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&trailer,sizeof(int),1);
//	    fprintf(stderr,"h %d t %d\n",header,trailer);
	    assert(header == trailer);
	    /*
	    ** hydro properties
	    */
	    fread(&header,sizeof(int),1,ad.GasFile);
	    if (ad.doswap) reorder(&header,sizeof(int),1);
	    assert(ad.Nhydroproperties == header/(sizeof(float)*ad.Ncell[i]));
	    for (j = 0; j < ad.Ncell[i]; j++) {
		assert(fread(cellhydroproperties,sizeof(float),ad.Nhydroproperties,ad.GasFile) == ad.Nhydroproperties);
		if (ad.doswap) reorder(cellhydroproperties,sizeof(float),ad.Nhydroproperties);
		}
	    assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
	    if (ad.doswap) reorder(&trailer,sizeof(int),1);
	    assert(header == trailer);
//	    fprintf(stderr,"h %d t %d\n",header,trailer);
	    /*
	    ** other properties
	    */
	    if (ad.GRAVITY || ad.RADIATIVE_TRANSFER) {
		assert(fread(&header,sizeof(int),1,ad.GasFile) == 1);
		if (ad.doswap) reorder(&header,sizeof(int),1);
		assert(ad.Notherproperties == header/(sizeof(float)*ad.Ncell[i]));
		for (j = 0; j < ad.Ncell[i]; j++) {
		    assert(fread(cellotherproperties,sizeof(float),ad.Notherproperties,ad.GasFile) == ad.Notherproperties);
		    if (ad.doswap) reorder(cellotherproperties,sizeof(float),ad.Notherproperties);
		    }
		assert(fread(&trailer,sizeof(int),1,ad.GasFile) == 1);
		if (ad.doswap) reorder(&trailer,sizeof(int),1);
		assert(header == trailer);
//		fprintf(stderr,"h %d t %d\n",header,trailer);
		}
	    }



/*
	for (i = 0; i < ad.Nhydroproperties; i++) {
	    fprintf(stderr,"i %d %g\n",i,cellhydroproperties[i]);
	    }
*/

	for (i = 0; i < ART_NUMBER_GAS_LEVELS; i++) {
	    fprintf(stderr,"i %d N %lu\n",i,ad.Ncell[i]);
	    }

	/*
	** Dark Matter and Stars
	*/
	ad.CoordinatesDataFile = fopen(ad.CoordinatesDataFileName,"r");
	assert(ad.CoordinatesDataFile != NULL);
	ac = realloc(ac,ad.Nparticleperrecord*sizeof(COORDINATES));
	assert(ac != NULL);
	pdp = realloc(pdp,gi.Nparticleperblockdark*sizeof(PROFILE_DARK_PARTICLE));
	assert(pdp != NULL);
	psp = realloc(psp,gi.Nparticleperblockstar*sizeof(PROFILE_STAR_PARTICLE));
	assert(psp != NULL);
	Nparticleread = 0;
	Nrecordread = 0;
	Icurrentblockdark = 0;
	Icurrentblockstar = 0;
	read_art_header_star(&ad);
	for (i = 0; i < ad.Nrecord; i++) {
	    read_art_coordinates_record(ad,ac);
	    fprintf(stderr,"i %d Nrecord %d\n",i,ad.Nrecord);
	    for (j = 0; j < ad.Nparticleperrecord; j++) {
		if (Nparticleread < ad.Ndark) {
		    /*
		    ** Dark Matter
		    */
		    for (k = 0; k < 3; k++) {
			pdp[Icurrentblockdark].r[k] = ac[j].r[k] + ad.shift;
			pdp[Icurrentblockdark].v[k] = ac[j].v[k];
			}
		    for (k = ad.Lmaxdark; k >=0; k--) {
			if (ad.ah.num[k] >= Nparticleread) L = ad.Lmaxdark-k;
			}
		    pdp[Icurrentblockdark].mass = ad.massdark[L];
		    Nparticleread++;
		    Icurrentblockdark++;
		    if ((Icurrentblockdark == gi.Nparticleperblockdark) || (Nparticleread == ad.Ndark)) {
			/*
			** Block is full or we reached end of dark matter particles
			*/
			gi.Nparticleinblockdark = Icurrentblockdark;
{'fraction_non_alphanumeric': 0.1768081694657777, 'fraction_numerical': 0.018246265507637777, 'mean_word_length': 2.8523961245854736, 'pattern_counts': {'":': 0, '<': 62, '<?xml version=': 0, '>': 443, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 4, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 1, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '438713', 'n_tokens_mistral': 24408, 'n_tokens_neox': 22879, 'n_words': 3985}