raybundle.c

用于2维的射线追踪

	bi->nb_bundle++;
	bi->bundle_set = (struct bundle_set_t *) 
		realloc (bi->bundle_set, sizeof(struct bundle_set_t)*bi->nb_bundle);
	assert(bi->bundle_set);

	/* filled */
	source.lat   = rd->src.lat;    source.lon = rd->src.lon;    
	source.prof = rd->src.prof;
	receiver.lat = rd->dest.lat; receiver.lon = rd->dest.lon; 
	receiver.prof = rd->dest.prof;

	bs = bi->bundle_set;
	bs[bi->nb_bundle - 1].receiver_lcid = lin_id_dest;
	bs[bi->nb_bundle - 1].bundle.source = dup_geo_coord(&source);
	bs[bi->nb_bundle - 1].bundle.receiver = dup_geo_coord(&receiver);
	bs[bi->nb_bundle - 1].bundle.ray_travel_time = rd->ray_travel_time;
	bs[bi->nb_bundle - 1].bundle.phase = strdup(rd->phase);
	bs[bi->nb_bundle - 1].nb_ray = 1;

	/* sort the  bundle set */
	qsort (&(bs[0]), bi->nb_bundle, sizeof(struct bundle_set_t), bundle_compare);
#ifdef DEBUG_RAYBUNDLE
	bundle_set_show(bi);
#endif	
}

void bundle_add_ray (struct bundle_info_t ****bi, 
		     struct raydata_t *raydata, 
		     struct mesh_t *mesh,
		     int ray_id) 
{
	   int res; 
	   int lin_id_dest; 
	   struct coord_z3_t cell_id_src; 
	   struct coord_z3_t cell_id_dest;
	   int x,y,z;

	   /* source */
    	   res = mesh_cellinfo_point2cell(mesh, 
			    raydata->src.lat, 
			    raydata->src.lon, 
			    raydata->src.prof, 
			    &cell_id_src);
	   if (res<=0) {
	   	fprintf(stderr, "ray(%d) source is not in the mesh\n", ray_id);
		print_raydata(stderr, raydata);
	   	exit (1);
	   }
	   x=cell_id_src.x;
	   y=cell_id_src.y;
	   z=cell_id_src.z;

	   /* dest */
	   res = mesh_cellinfo_point2cell(mesh, 
			    raydata->dest.lat, 
			    raydata->dest.lon, 
			    raydata->dest.prof, 
			    &cell_id_dest);
	   if (res<=0) {
	   	fprintf(stderr, "ray(%d) dest is not in the mesh\n", ray_id);
		print_raydata(stderr, raydata);
	   	exit (1);
	   }
    	   lin_id_dest = linearize_cell_id (&cell_id_dest, mesh);

	   /* alloc if needed */
	   if (!bi[z][x][y]) {
		   bi[z][x][y] = (struct bundle_info_t *)
			   malloc (sizeof(struct bundle_info_t));
		   assert(bi[z][x][y]);
		   bi[z][x][y]->bundle_set = NULL;
		   bi[z][x][y]->nb_bundle = 0;
	   }

	   /* feed bundle info */
#ifdef DEBUG_RAYBUNDLE
	   fprintf(stderr, "bundle_add_ray: bi[%d][%d][%d], lin_id_dest=%d, phase='%s'\n", 
			   z,x,y,lin_id_dest,raydata->phase);
#endif
	   bundle_set_insert (bi[z][x][y], lin_id_dest, raydata);
	   
}

/** \brief Write ray bundle to file 
 */
void bundle_write (struct bundle_info_t ****bi, 
		   struct mesh_t *m,
		   FILE *bundle_fd) 
{
    int nb_lat_cell, nb_lon_cell;
    int z, x, y, n;
    struct bundle_t *bundle;

    for (z = 0; z < m->nlayers; z++) {
	nb_lat_cell = m->layer[z]->nlat;
	nb_lon_cell = m->layer[z]->nlon;

        for (x = 0; x < nb_lat_cell; x++) { 
	    for (y = 0; y < nb_lon_cell; y++) { 	
		if (!bi[z][x][y]) {
			continue;
		}
		for (n=0; n < bi[z][x][y]->nb_bundle; n++) {
		    bundle = &(bi[z][x][y]->bundle_set[n].bundle);

	            fprintf(bundle_fd, "%d %f %f %f %f %f %f %s %f\n",
				    bi[z][x][y]->bundle_set[n].nb_ray,
				    bundle->source->lat, 
				    bundle->source->lon, 
				    bundle->source->prof,
				    bundle->receiver->lat, 
				    bundle->receiver->lon, 
				    bundle->receiver->prof,
				    bundle->phase,
				    bundle->ray_travel_time); 
		}
	    }
	}
    }

}

/** \brief write into sco file the number of ray bundle per cell
 *
 * We are exporting data to sco format from a light mesh.
 * The next step is probably to re-import this sco data file into
 * a real mesh to obtain a vtk file. 
 *
 */
void bundle_write_sco (struct bundle_info_t ****bi, 
		   struct mesh_t *m,
		   FILE *sco_fd) 
{
    int nb_lat_cell, nb_lon_cell;
    int z, x, y;
    int nb_ray, n;

    fprintf(sco_fd, "# format=sco, generated by %s v%s, mesh=%s\n", 
		    RAYBUNDLE_PACKAGE, RAYBUNDLE_VERSION, m->xml_filename);
    fprintf(sco_fd, "%d nb_of_bundles nb_of_rays\n", 2);

    for (z = 0; z < m->nlayers; z++) {
        nb_lat_cell = m->layer[z]->nlat;
        nb_lon_cell = m->layer[z]->nlon;

        for (x = 0; x < nb_lat_cell; x++) { 
	    for (y = 0; y < nb_lon_cell; y++) { 	
		if (!bi[z][x][y]) {
			continue;
		}

		nb_ray = 0;
		for (n=0; n < bi[z][x][y]->nb_bundle; n++) {
			nb_ray += bi[z][x][y]->bundle_set[n].nb_ray;
		}

		/* coordinates in the real mesh */
	        fprintf(sco_fd, "[%d,%d,%d] %d %d\n", 
			     x, y, z, bi[z][x][y]->nb_bundle, nb_ray);
	    }
	}
    }

}







int main(int argc, char **argv)
{ 
    char *mesh_file = NULL;
    char *ray_file = NULL;
    char *bundle_file = NULL;
    char *sco_file = NULL;

    struct mesh_t *mesh;
    FILE *ray_fd;
    FILE *bundle_fd;
    FILE *sco_fd=NULL;

    int nbread = 0;
    int nberr = 0;
    int SIZE = 1000;
    int ray_id = 0;

    int i;
    struct raydata_t *raydata;
    struct bundle_info_t ****bundle_info = NULL;
    
    /* progress bar */
    int max;
    int cumulated_size = 0;
    


    parse_command_line(argc, argv, 
		       &mesh_file, &ray_file, &bundle_file, &sco_file);
	
    /* init of the mesh->parameters struct */
    if (!(mesh = mesh_init_from_file(mesh_file))) {
	exit(1);
    }

    /* check if sco section already available */
    if (sco_file && mesh->data[SCO] && mesh->data[SCO]->ndatafile) {
	    fprintf(stderr, "Detected an already sco file in the mesh '%s'\n",
			    mesh_file);
	    exit(1);
    }

#ifdef DEBUG_RAYBUNDLE
    mesh_dump_parameter(stdout, "", "", mesh->parameter);
#endif

    /* open rays input file */
    if (!(ray_fd = fopen(ray_file, "r"))) {
	fprintf(stderr, "Can not open input ray file '%s' ... exiting!\n",
			ray_file);
	exit(1);
    }

    /* open raybundle output file */
    if (!(bundle_fd = fopen(bundle_file, "w"))) {
	fprintf(stderr, "Can not open bundle ray file '%s' ... exiting!\n",
			bundle_file);
	exit(1);
    }

    /* open sco output file */
    if (sco_file) {
        if (!(sco_fd = fopen(sco_file, "w"))) {
	    fprintf(stderr, "Can not open sco file '%s' ... exiting!\n",
			sco_file);
	    exit(1);
        }
    }

    
    /* create the structure  */
    bundle_info = bundle_alloc(mesh);
    assert(bundle_info);

    fprintf(stdout, "counting number of ray ... ");
    fflush(stdout);
    max = get_number_of_lines(ray_fd);
    fprintf(stdout, "%d rays\n", max);

    fprintf(stdout, "making bundle ... ");
#ifndef DEBUG_RAYBUNDLE
    fflush(stdout);
    start_progress_bar_percent();
#else
	fprintf(stdout, "\n");
#endif

    while (!feof(ray_fd)) {
#ifndef DEBUG_RAYBUNDLE
	progress_bar_percent(cumulated_size, max);
#endif	
    	raydata = get_raydata(ray_fd, SIZE, &nbread, &nberr);
	cumulated_size+=SIZE;
	
	for (i = 0; i < nbread; i++) {
	    /* raydata is in radian, convert it to DEG */
    	    raydata[i].src.lat *= TO_DEG;
	    raydata[i].src.lon *= TO_DEG;
	    raydata[i].dest.lat *= TO_DEG;
	    raydata[i].dest.lon *= TO_DEG;
		
	   /* take care of longitude sign, mesh lon_min is defined > 0 */
    	   if ( raydata[i].src.lon < 0. || mesh->parameter->lon_max > 360.) {
		   raydata[i].src.lon += 360.;
	   }
	   if ( raydata[i].dest.lon < 0. || mesh->parameter->lon_max > 360.) {
		   raydata[i].dest.lon += 360.;
	   }

	   /* update bundle info */
	   ray_id++;
	   bundle_add_ray (bundle_info, &(raydata[i]), mesh, ray_id); 
	}
    }
    fclose (ray_fd);
    
#ifndef DEBUG_RAYBUNDLE
    stop_progress_bar_percent();
#endif

    /* write results and some stats */
    fprintf(stdout, "writing bundle file '%s'\n", bundle_file);
    bundle_write(bundle_info, mesh, bundle_fd);
    fclose (bundle_fd);
    
    /* generate sco data file */
    if (sco_file) {
	fprintf(stdout, "writing sco file '%s'\n", sco_file);
        bundle_write_sco (bundle_info, mesh, sco_fd);
        mesh_add_data_filename(mesh, SCO, sco_file);
        fclose (sco_fd);
	
        /* update the mesh xml file with the new sco section */
    	mesh2xml(mesh, mesh_file);
    }

    return(0);
}