raybundle.c

用于2维的射线追踪

#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);
}

void
parse_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 */