/* * Ray2mesh : software for geophysicists. * Compute various scores attached to the mesh cells, based on geometric   information that rays bring when the traverse the cell. * * Copyright (C) 2003, St閜hane Genaud and Marc Grunberg * * This tool is free software; you can redistribute it and/or * modify it under the terms of the GNU Library General Public * License as published by the Free Software Foundation; either * version 2 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU * Library General Public License for more details. * * You should have received a copy of the GNU Library General Public * License along with this library; if not, write to the Free * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */#ifdef HAVE_CONFIG_H#include <config.h>#endif#ifdef USE_MPI#include <mpi.h>#endif#ifdef USE_ZLIB#include "zlib.h"#endif#include <const_ray2mesh.h>#include <stdio.h>#include <stdlib.h>#include <string.h>#include <unistd.h>#include <sys/time.h>#include <sys/types.h>#include <sys/stat.h>#include <errno.h>#include <signal.h>#include <popt.h>#include <iasp/iasp.h>#include <ray/raydescartes.h>#include <mesh/mesh.h>#include <mesh/layer.h>#include <mesh/cellinfo.h>#include <mesh/mesh2xml.h>#include <sparse/sparse.h>#include "ray2mesh.h"#include "r2mfile.h"#include "merge.h"#include "score.h"#include "buffer.h"#include "filter.h"char *celldatafilename = NULL;		/**< these are kept global to close */FILE *fdinput = NULL;			/**< them in case of emergency (Ctrl-C) */struct mesh_t *mesh;			/**< mesh structure */struct cell_info_t ****cell_info = NULL;/**< traversed cells information */char *output_format = NULL;			/**< r2m and/or sco output */int rank;volatile int mem_used=0;/** \brief shows ray2mesh version ans compilation flags */void ray2mesh_info(struct ray_config_t *ray_config){    char *modelver = libmodelversion();    char *rayver = librayversion();    char *raycodever = libraycodeversion();    char *meshver = libmeshversion();    char *sparsever = libsparseversion();    fprintf(stderr, "%s: version %s", PACKAGE, VERSION);#ifdef STATIC_BIN    fprintf(stderr, " [*static binary*]");#endif    fprintf(stderr, "\nusing \t%s\n", modelver);    fprintf(stderr, "\t%s\n", rayver);    fprintf(stderr, "\t%s\n", raycodever);    fprintf(stderr, "\t%s\n", meshver);    fprintf(stderr, "\t%s\n", sparsever);    free(modelver);    free(rayver);    free(raycodever);    free(meshver);    free(sparsever);    fprintf(stderr, "Compilation flags : ");#ifdef USE_MPI    fprintf(stderr, "USE_MPI ");    #ifdef MASTER_SLAVE        fprintf(stderr, "MASTER_SLAVE ");    #else        fprintf(stderr, "SCATTER ");    #endif    #ifdef USE_MPILB    fprintf(stderr, "USE_MPILB ");    #endif#endif#ifdef DEBUG    fprintf(stderr, "DEBUG ");#endif#ifdef ALLOC_MEMORY_CHUNK    fprintf(stderr, "ALLOC_MEMORY_CHUNK ");#endif#ifdef RAYTRACING_ONLY    fprintf(stderr, "RAYTRACING_ONLY ");#endif#ifdef USE_ZLIB    fprintf(stderr, "USE_ZLIB(%s) ", zlibVersion());#endif    fprintf(stderr, "\n");}/** \brief given a filename, stops if file is unreadable  */void check_file_access(char *filename){struct stat st;    /* check read permissions */    if (access(filename, R_OK) < 0) {	fprintf(stderr, "%s : cannot access file %s. Exiting.\n", PACKAGE, filename);	exit(1);    }    /* check filetype */    stat(filename,&st);    if (!S_ISREG(st.st_mode)) {	fprintf(stderr, "%s : file %s is not a regular file. Exiting.\n", PACKAGE, filename);	exit(1);    }}/** \brief handler for Ctrl-C, close the files */void emergency_halt(){    fprintf(stdout, "*** Emergency stop (Ctrl-C). Sync files ***\n");    fflush(stdout);    if (strchr(output_format, 'r')) {	/* R2M */	make_domain_info_file("emergency.r2m", cell_info, mesh, 0, 1, 0);    } else {	/* SCO */        fprintf(stdout, "-> Computing cell score ... ");	compute_score(cell_info, mesh);        fprintf(stdout, "done.\n-> Dumping cell information ... ");	mesh_cellinfo_write_sco ("emergency.sco", cell_info, mesh);    }    fprintf(stdout, "done.\n");#ifdef USE_MPI    MPI_Finalize();#endif    exit(1);}/* * \brief parse command line using popt : checks validity of optoins  * specified and returns parameters values. * * PRECOND : the default values for optional arguments must be set.      */voidparse_command_line(int argc, char **argv,		   struct ray_config_t *ray_config,		   struct ray_filter_t *ray_filter,		   char **meshfile,		   char **outputfile,		   char **inputdatafile,		   char **rayfilter_filename, 		   float *mem_limit, 		   char **output_format,		   char **tmpdir#ifdef MASTER_SLAVE		   , int *chunksize#endif    ){    static float length_opt;    static char *meshfile_opt=NULL;	/* there to build on IRIX with popt */    static char *outputfile_opt=NULL;    static char *vmodel;    static char *inputdatafile_opt=NULL;    static char *iasptables_opt;    static char *filtered=NULL;    static float mem_limit_opt;    static char *setfilter_opt=NULL;    static char *outputformat_opt=NULL;    static int iterate_opt;    static char *tmpdir_opt=NULL;#ifdef MASTER_SLAVE    static int chunksize_opt;#endif        poptContext cx;    int rc;    int nbarg;    #include "options.h"    /* init */    set_ray_config_to_default(ray_config);    /* ray_config->iterative_mode = ITERATIVE_MODE_ANGULAR; */    /* ray_config->iterative_mode = ITERATIVE_MODE_OFF; */    /* ray_config->iterative_mode = ITERATIVE_MODE_ON; */    /* args parsing using popt */    cx = poptGetContext(PACKAGE, argc, (const char **) argv, options, 0);    if (argc == 1) {	poptPrintUsage(cx, stdout, 0);	exit(1);    }    do {	rc = poptGetNextOpt(cx);	switch (rc) {	    /* help */	case OPT_HELP:	    poptPrintHelp(cx, stdout, 0);	    exit(0);	case OPT_VELOCITY_MODEL:            check_file_access(vmodel);	    ray_config->velocity_model = load_velocity_model(vmodel);	    ray_config->get_velocity_p = my_get_velocity_p;	    ray_config->get_velocity_s = my_get_velocity_s;	    break;	case OPT_LENGTH:	    ray_config->length_inc = length_opt;	/* update if increment specified */	    break;	case OPT_OUTPUT_FORMAT:	/* sco output selected */	    if ( !strchr(outputformat_opt,'s') &&                 !strchr(outputformat_opt,'r'))            {                    fprintf(stderr, "format must be 'sco' or 'r2m'.\n");                    exit(1);            }	    if (strchr(outputformat_opt,'s') && strchr(outputformat_opt,'r')) {	            fprintf(stderr, "sorry, both output not yet implemented\n");                    exit(1);	    }	    break;#ifdef MASTER_SLAVE	case OPT_CHUNKSIZE:	    *chunksize = chunksize_opt;	    break;#endif	case OPT_TABLES:#ifdef DEBUG	    fprintf(stderr, "setting iasp tables path to \"%s\"\n",		    iasptables_opt);#endif	    set_modnam(iasptables_opt);	    break;	case OPT_ITERATE:	    ray_config->iterative_mode = ITERATIVE_MODE_ON;	    break;	case OPT_MEM_LIMIT:	    *mem_limit = mem_limit_opt;	    break;	case OPT_SET_FILTER:	    nbarg = parse_filter_value(strdup(setfilter_opt), ray_filter);	    ray_filter->activated = 1;	    if (nbarg!=3) {		    fprintf(stderr, "parsing filter string (%s) failed\n", 				    setfilter_opt); 		    fprintf(stderr, "use --set-filter=r,d,n\n");		    exit(1);	    }	    break;	case OPT_VERSION:	    ray2mesh_info(ray_config);	    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);    /* check meshfile */    if (!meshfile_opt) {	fprintf(stderr, "%s: missing -m option\n", PACKAGE);	exit(1);    } else {	*meshfile = strdup (meshfile_opt);    }    if (!outputfile_opt) {        fprintf(stderr, "%s: missing -o option\n", PACKAGE);	exit(1);    } else {    	*outputfile = strdup(outputfile_opt);    }    /* check inputdatafile option */    if (!inputdatafile_opt) {	fprintf(stderr, "%s: missing -i option\n", PACKAGE);	exit(1);    } else {    	*inputdatafile = strdup(inputdatafile_opt);    }    /* check output format */    if (!outputformat_opt) {        fprintf(stderr, "%s: missing -f option, sco or r2m\n", PACKAGE);	exit(1);    } else {    	*output_format = strdup(outputformat_opt);    }    if (filtered) {    	*rayfilter_filename = strdup (filtered);    }    /* set tmp directory */    if (tmpdir_opt) {    	*tmpdir = strdup (tmpdir_opt);    } else {	*tmpdir = strdup ("/tmp");    }	    {     struct stat stabuf;     if (stat(*tmpdir, &stabuf) == -1) {        perror(*tmpdir);	exit(1);     }     if (!S_ISDIR(stabuf.st_mode)) {	fprintf(stderr, "%s is not a directory !\n", *tmpdir);	exit(1);     }    }    poptFreeContext(cx);}void print_ray_ko_info (int rank, struct raydata_t *raydata, int i, int nbrays, char *rai_status) {	fprintf(stderr,		   "[process %d] KO %d/%d,  %s,  -s %f,%f,%f -d %f,%f,%f -p %s\n",		   rank,		   i,nbrays,		   rai_status?rai_status:"",		   raydata[i].src.lat * TO_DEG,		   raydata[i].src.lon * TO_DEG,		   raydata[i].src.prof,		   raydata[i].dest.lat * TO_DEG,		   raydata[i].dest.lon * TO_DEG,		   raydata[i].dest.prof, raydata[i].phase); }/* \brief given a set of raydata elements (src and dest) pairs  * for rays, and the number of rays contained in the set,                * run ray3D_descartes of this bunch of ray data with a discretization   * step length of length_inc km. nb_ray_computed and nb_ray_rejected are * updated accordingly                                                   *  * POSTCOND : raydata is freed before returning                          */void bunch_of_ray(struct raydata_t *raydata,		  const int nbrays,		  const long int offset,		  struct ray_config_t *ray_config,		  struct ray_filter_t *ray_filter,		  int *nb_ray_computed,		  int *nb_ray_rejected,		  FILE * filterfd,		  FILE * matrix_length_fd, 		  FILE * residual_time_fd,		  FILE * event_fd){    struct ray3D_t **my3Drays, **ptr;    int i;    char **phases, **cp_phases;    char *rai_status=NULL;    double dT;			/* time residual */    double delta;    fprintf(stdout, "*** [process %d] bunch_of_ray rays  rayId=[%d  --> %d]\n",	    rank, *nb_ray_computed, *nb_ray_computed + nbrays -1);    for (i = 0; i < nbrays; i++) {	if (ray_filter->activated) {	        /* check if nb bundle is ok */		if (raydata[i].event_id < ray_filter->nb_event_min) {			print_ray_ko_info (rank, raydata, i, nbrays, "nbbundle");			(*nb_ray_rejected)++;			continue;		}	        /* check if delta is ok */		delta = AOB(&raydata[i].src, &raydata[i].dest)*TO_DEG;		if (delta > ray_filter->delta_max) {			print_ray_ko_info (rank, raydata, i, nbrays, "delta");			(*nb_ray_rejected)++;			continue;		}	}	/* do the ray-tracing */	phases = cp_phases = parse_phase(raydata[i].phase);	my3Drays = ray3D_descartes(phases, &raydata[i].src, 			&raydata[i].dest, ray_config);	/* desalloc phases */	while (*cp_phases) {	    free(*cp_phases);	    cp_phases++;	}	free(phases);	/* loop over all returned rays phases for this ray */	for (ptr = my3Drays; *ptr; ptr++) {	    /* dT = Tobserved - Tcomputed */	    dT = raydata[i].ray_travel_time - (*ptr)->time;	    /* check if residual time is small enough */	    if (ray_filter->activated && 	        fabs(dT) > ray_filter->residual_max) 	    {	    	/* free */		print_ray_ko_info (rank, raydata, i, nbrays, "residual");            	(*nb_ray_rejected)++;            	free_ray3D_descartes(*ptr);            	*ptr = NULL;		continue;	    }		    rai_status = str_ray_status((*ptr)->status);	    /* ray-trace failed */	    if ((*ptr)->status != RAYTRACE_OK) {		print_ray_ko_info (rank, raydata, i, nbrays, rai_status);		(*nb_ray_rejected)++;		free_ray3D_descartes(*ptr);            	*ptr = NULL;		free(rai_status);		rai_status=NULL;		continue;	    }	    /* ray-trace is ok */	    fprintf(stderr, "[process %d] OK ray %d/%d (phase %s) %s length=%f nbiter=%d dist_err=%f\n",		    rank,		    i, nbrays, raydata[i].phase, rai_status, 		    (*ptr)->total_length,		    (*ptr)->nb_iter_angle+(*ptr)->nb_iter_length,		    (*ptr)->dist_err);	    free(rai_status);	    rai_status=NULL;	    /* ray is ok */	    if (filterfd) {		/* user want to keep info about good rays */		fprintf(filterfd, "%ld %g %g %g %g %g %g %s %g %f\n",			raydata[i].event_id,			raydata[i].src.lat * TO_DEG,			raydata[i].src.lon * TO_DEG,			raydata[i].src.prof,			raydata[i].dest.lat * TO_DEG,			raydata[i].dest.lon * TO_DEG, 			raydata[i].dest.prof,			raydata[i].phase,			/*(*ptr)->code, */			raydata[i].ray_travel_time, 			dT);	    }#ifndef RAYTRACING_ONLY	    /* we must compute the time residual */	    if (residual_time_fd) {		fprintf(residual_time_fd, "%ld %lf\n",			(long int) *nb_ray_computed, (double) dT);	    }	    /* write event_id ray_id travel_time */   	    if (event_fd) {		fprintf(event_fd, "%ld %ld %lf\n", 			raydata[i].event_id, (long int) *nb_ray_computed, 			(*ptr)->time);	    }	    /* pass previous cell_info to accumulate results 	     * from several rays */	    cell_info = ray2mesh2(*ptr, *nb_ray_computed, mesh, cell_info,			 matrix_length_fd);#endif	    /* free */	    free_ray3D_descartes(*ptr);	    *ptr = NULL;	    (*nb_ray_computed)++;	}	free(my3Drays);    }    free(raydata);    }/*-----------------------------------------*//* main                                    *//*------- ---------------------------------*/#ifdef USE_MPI#ifdef MASTER_SLAVE#include "main_mpi_ms.c"#else#include "main_mpi.c"#endif#else#ifdef DEBUG#include <mcheck.h>#endif#include "main_seq.c"#endif