#include <stdio.h>#include <unistd.h>#include <string.h>#include <stdlib.h>#include <popt.h>#include <mesh/mesh.h>#include <mesh/layer.h>#include <mesh/cell.h>#include <mesh/cellinfo.h>#include <mesh/mesh2xml.h>#include <raycode/raycode.h>#include <ray/raydescartes.h>#include "buffer.h"#include "ray2mesh.h"#define RAYBUNDLE_VERSION "0.2"#define RAYBUNDLE_PACKAGE "raybundle"/*#define DEBUG_RAYBUNDLE*/struct bundle_t {	struct coord_geo_t *source;	struct coord_geo_t *receiver;	double ray_travel_time;	char *phase;};struct bundle_set_t {	int receiver_lcid;	int nb_ray;	struct bundle_t bundle;};struct bundle_info_t {	int nb_bundle;	struct bundle_set_t * bundle_set;};/* --- duplicated progress bar from gimbos --- */void start_progress_bar_percent (void) {        fprintf(stdout, "      ");        fflush(stdout);}float progress_bar_percent (int nb, int max) {        float percent=100.*(float)nb/(float)max;        fprintf(stdout, "\b\b\b\b\b\b%5.1f%%", percent);        fflush(stdout);        return(percent);}void stop_progress_bar_percent (void) {        fprintf(stdout, "\b\b\b\b\b\bcompleted\n");        fflush(stdout);}/* ------------------------------------------- */void raybundle_version (void) {    char *meshver = libmeshversion();        fprintf(stderr, "raybundle: version %s", RAYBUNDLE_VERSION);    fprintf(stderr, "\nusing \t%s\n", meshver);    free(meshver);}voidparse_command_line(int argc, char **argv,                   char **mesh_file,                   char **ray_file,                   char **bundle_file,		   char **sco_file){    static char *mesh_file_opt = NULL;  /* there to build on IRIX with popt */    static char *ray_file_opt = NULL;    static char *bundle_file_opt = NULL;    static char *sco_file_opt = NULL;    poptContext cx;    int rc;    #include "raybundle_options.h"	    /* args parsing using popt */    cx = poptGetContext(RAYBUNDLE_PACKAGE, argc, 		       (const char **) argv, options, 0);    if (argc == 1) {        poptPrintUsage(cx, stdout, 0);        exit(1);    }    do {        rc = poptGetNextOpt(cx);        switch (rc) {        case OPT_HELP:		/* help */            poptPrintHelp(cx, stdout, 0);            exit(0);	case OPT_BUNDLE_FILE:	    if (!access(bundle_file_opt, F_OK)) {	        fprintf(stderr, "file '%s' already exists !\n", bundle_file_opt);		exit(1);	    }	    break;	case OPT_SCO_FILE:	    if (!access(sco_file_opt, F_OK)) {	        fprintf(stderr, "file '%s' already exists !\n", sco_file_opt);		exit(1);	    }	    break;	case OPT_VERSION:	    raybundle_version();	    exit(0);	case POPT_ERROR_BADOPT: /* error */            fprintf(stderr, "%s: %s\n", poptBadOption(cx, 0),                    poptStrerror(rc));            poptPrintUsage(cx, stdout, 0);            exit(1);        }    }    while (rc > 0);    /* bundle file */    *bundle_file = bundle_file_opt;    if (!*bundle_file) {        fprintf(stderr, "%s: missing -o option\n", RAYBUNDLE_PACKAGE);        exit(1);    }    /* sco file */    *sco_file = sco_file_opt;   /* check meshfile */    *mesh_file = mesh_file_opt;     if (!*mesh_file) {        fprintf(stderr, "%s: missing -m option\n", RAYBUNDLE_PACKAGE);        exit(1);    }    /* check inputdatafile option */    *ray_file = ray_file_opt;    if (!*ray_file) {        fprintf(stderr, "%s: missing -i option\n", RAYBUNDLE_PACKAGE);        exit(1);    }    poptFreeContext(cx);}void print_raydata (FILE *fd, struct raydata_t *rd) {	fprintf(fd, "%f %f %f %f %f %f %s\n",			rd->src.lat, rd->src.lon, rd->src.prof, 			rd->dest.lat, rd->dest.lon, rd->dest.prof, 			rd->phase);}/** \brief Allocation for the information bundle structure  **/struct bundle_info_t ****bundle_alloc(struct mesh_t *m) {	struct bundle_info_t ****bi;	int z, x, y;    	int nb_lat_cell, nb_lon_cell;	bi = (struct bundle_info_t ****) 		malloc(m->nlayers*sizeof(struct bundle_info_t ***));    	assert(bi);    	for (z = 0; z < m->nlayers; z++) {		nb_lat_cell = m->layer[z]->nlat;		nb_lon_cell = m->layer[z]->nlon;		bi[z] = (struct bundle_info_t ***) 			malloc(nb_lat_cell * sizeof(struct bundle_info_t **));		assert(bi[z]);		for (x = 0; x < nb_lat_cell; x++) { 			bi[z][x] = (struct bundle_info_t **) 				malloc(nb_lon_cell * sizeof(struct bundle_info_t *)); 			assert(bi[z][x]); 						for (y = 0; y < nb_lon_cell; y++) { 				bi[z][x][y] = NULL; 			} 		}	}	return (bi);}int bundle_compare (const void * a, const void * b){    const struct bundle_set_t *aa = (const struct bundle_set_t *) a;    const struct bundle_set_t *bb = (const struct bundle_set_t *) b;    int val;    if ((*aa).receiver_lcid  > (*bb).receiver_lcid) {       return 1;    } else if ( (*aa).receiver_lcid  == (*bb).receiver_lcid) {	   val = strcmp((*aa).bundle.phase,(*bb).bundle.phase);	   assert(val);	   return(val);    }    return -1;}int bundle_update (struct bundle_set_t *bs, struct raydata_t *rd, char *mean_phase) {	int nb = bs->nb_ray;		/* update source */	bs->bundle.source->lat = 		(nb*bs->bundle.source->lat + rd->src.lat)/(nb+1);	bs->bundle.source->lon = 		(nb*bs->bundle.source->lon + rd->src.lon)/(nb+1);	bs->bundle.source->prof = 		(nb*bs->bundle.source->prof + rd->src.prof)/(nb+1);		/* update receiver */	bs->bundle.receiver->lat = 		(nb*bs->bundle.receiver->lat + rd->dest.lat)/(nb+1);	bs->bundle.receiver->lon = 		(nb*bs->bundle.receiver->lon + rd->dest.lon)/(nb+1);	bs->bundle.receiver->prof = 		(nb*bs->bundle.receiver->prof + rd->dest.prof)/(nb+1);	/* update travel time */	bs->bundle.ray_travel_time = 		(nb*bs->bundle.ray_travel_time + rd->ray_travel_time)/(nb+1);	/* update to the mean phase */	free(bs->bundle.phase);	bs->bundle.phase = strdup (mean_phase);		bs->nb_ray++;	return(bs->nb_ray);}void bundle_set_show (struct bundle_info_t *bi) {	int i;		for (i=0; i<bi->nb_bundle; i++) {		fprintf(stderr, "    bundle[%d] : lincelid=%d ... nbrays=%d\n", 				i, bi->bundle_set[i].receiver_lcid, 				bi->bundle_set[i].nb_ray);	}}/** \brief Insert a ray into a bundle * * If the current ray is part of a pre-existing bundle,  * the corresponding bundle is updated accordingly. * If not a new bundle is created and initialized. * The bundle_set array is sorted by the linear destination cell id as first * key, and phase as second key. *  */void bundle_set_insert (struct bundle_info_t *bi, 			int lin_id_dest,			struct raydata_t *rd){		int i=0;	struct bundle_set_t *bs;	struct coord_geo_t source;	struct coord_geo_t receiver;	char *ray_mean;		if (bi->nb_bundle) {		/* find a previous existing bundle */		bs = bi->bundle_set;		for (i=0; i<bi->nb_bundle; i++) {			if (lin_id_dest < bs[i].receiver_lcid) {				break;			}			if ( lin_id_dest == bs[i].receiver_lcid  &&			     (ray_mean = raycode_compute_mean (rd->phase, 1, bs[i].bundle.phase, bs[i].nb_ray)) != NULL )			{				int nb;				/*fprintf(stderr, "ray_mean='%s'\n", ray_mean);*/				/* update this bundle */				nb = bundle_update( &(bs[i]), rd, ray_mean);				free(ray_mean);#ifdef DEBUG_RAYBUNDLE				fprintf(stderr, "  bundle_set_insert: found previous bundle, %d rays\n",						nb);				bundle_set_show(bi);#endif				return;			}		}	}#ifdef DEBUG_RAYBUNDLE	fprintf(stderr, "  bundle_set_insert: no previous bundle\n");#endif		/* add a new bundle at the end of bundle set */