Changeset 5485


Ignore:
Timestamp:
Jul 10, 2008, 2:14:22 PM (11 years ago)
Author:
ole
Message:

Reverted MUX reader to revision 5470 - the last working version.
This one still has servere memory leaks and needs to be redeveloped.

Location:
anuga_core/source/anuga/shallow_water
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r5474 r5485  
    4848
    4949"""
     50
     51# This file was reverted from changeset:5484 to changeset:5470 on 10th July
     52# by Ole.
    5053
    5154import exceptions
     
    48384841NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'   
    48394842
    4840 def read_mux2_py(filenames,weights,permutation=None,verbose=False):
    4841     """
    4842     Access the mux files specified in the filenames list. Combine the
    4843     data found therin as a weighted linear sum as specifed by the weights.
    4844     If permutation is None extract timeseries data for all gauges within the
    4845     files.
    4846    
    4847     Input:
    4848     filenames:   List of filenames specifiying the file containing the
    4849                  timeseries data (mux2 format) for each source
    4850     weights:     Weights associated with each source
    4851     permutation: Specifies the gauge numbers that for which data is to be
    4852                  extracted
    4853     """
     4843def read_mux2_py(filenames,weights,verbose=False):
     4844
    48544845    from Numeric import ones,Float,compress,zeros,arange
    48554846    from urs_ext import read_mux2
     
    48654856        verbose=0
    48664857       
    4867     data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
     4858    data=read_mux2(numSrc,filenames,weights,file_params,verbose)
    48684859
    48694860    msg='File parameter values were not read in correctly from c file'
     
    50465037                raise IOError, msg
    50475038
    5048     if ordering_filename is not None:
    5049         # Read ordering file
    5050        
    5051         try:
    5052             fid=open(ordering_filename, 'r')
    5053             file_header=fid.readline().split(',')
    5054             ordering_lines=fid.readlines()
    5055             fid.close()
    5056         except:
    5057             msg = 'Cannot open %s'%ordering_filename
    5058             raise Exception, msg
    5059 
    5060         reference_header = 'index, longitude, latitude\n'
    5061         reference_header_split = reference_header.split(',')
    5062         for i in range(3):
    5063             if not file_header[i].strip() == reference_header_split[i].strip():
    5064                 msg = 'File must contain header: '+reference_header+'\n'
    5065                 raise Exception, msg
    5066 
    5067         if len(ordering_lines)<2:
    5068             msg = 'File must contain at least two points'
    5069             raise Exception, msg
    5070 
    5071    
    5072        
    5073         permutation = [int(line.split(',')[0]) for line in ordering_lines]
    5074     else:
    5075         permutation = None 
    5076 
     5039    #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
     5040    #for now set x_momentum and y_moentum quantities to zero
    50775041    if (verbose): print 'reading mux2 file'
    50785042    mux={}
    50795043    for i,quantity in enumerate(quantities):
    50805044        # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
    5081         times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
     5045        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
    50825046        if quantity!=quantities[0]:
    50835047            msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
     
    50925056        elevation_old=elevation
    50935057        starttime_old=starttime
     5058
    50945059       
    50955060    if ordering_filename is not None:
    5096         latitudes = latitudes_urs
    5097         longitudes = longitudes_urs
     5061        # Read ordering file
     5062       
     5063        try:
     5064            fid=open(ordering_filename, 'r')
     5065            file_header=fid.readline().split(',')
     5066            ordering_lines=fid.readlines()
     5067            fid.close()
     5068        except:
     5069            msg = 'Cannot open %s'%ordering_filename
     5070            raise Exception, msg
     5071
     5072        reference_header = 'index, longitude, latitude\n'
     5073        reference_header_split = reference_header.split(',')
     5074        for i in range(3):
     5075            if not file_header[i] == reference_header_split[i]:
     5076                msg = 'File must contain header: '+reference_header+'\n'
     5077                raise Exception, msg
     5078
     5079        if len(ordering_lines)<2:
     5080            msg = 'File must contain at least two points'
     5081            raise Exception, msg
     5082
     5083   
     5084       
     5085        permutation = [int(line.split(',')[0]) for line in ordering_lines]
     5086   
     5087        latitudes = take(latitudes_urs, permutation)
     5088        longitudes = take(longitudes_urs, permutation)
    50985089       
    50995090        # Self check - can be removed to improve speed
     
    51285119        longitudes=longitudes_urs
    51295120        permutation = range(latitudes.shape[0])               
     5121       
    51305122
    51315123    msg='File is empty and or clipped region not in file region'
     
    51895181    if verbose: print 'Converting quantities'
    51905182    for j in range(len(times)):
    5191         for i in range(number_of_points):
    5192             w = zscale*mux['HA'][i,j] + mean_stage
    5193             h=w-elevation[i]
     5183        for i, index in enumerate(permutation):
     5184            w = zscale*mux['HA'][index,j] + mean_stage
     5185            h=w-elevation[index]
    51945186            stage[j,i] = w
    51955187
    5196             xmomentum[j,i] = mux['UA'][i,j]*h
    5197             ymomentum[j,i] = mux['VA'][i,j]*h
     5188            xmomentum[j,i] = mux['UA'][index,j]*h
     5189            ymomentum[j,i] = mux['VA'][index,j]*h
    51985190
    51995191    outfile.close()
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5473 r5485  
    22#
    33
     4# This file was reverted from changeset:5484 to changeset:5470 on 10th July
     5# by Ole.
    46
    57import unittest
     
    61056107        return base_name, files
    61066108
    6107     def test_read_mux2_py0(self):
     6109    def test_read_mux2_pyI(self):
    61086110        """Constant stage,momentum at each gauge
    61096111        """
     
    65436545        """Test multiple sources with ordering file
    65446546        """
     6547       
    65456548        tide=0
    65466549        time_step_count = 3
     
    69446947            fid.write(line)
    69456948        fid.close()
    6946        
     6949
    69476950        sts_file=base_name
    69486951        urs2sts(base_name, basename_out=sts_file,
     
    88888891if __name__ == "__main__":
    88898892
    8890     #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux2_py0')
    88918893    suite = unittest.makeSuite(Test_Data_Manager,'test')
    88928894    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5482 r5485  
    33gcc -shared urs_ext.o  -o urs_ext.so
    44*/
     5
     6/*
     7This file was reverted from changeset:5484 to changeset:5470 on 10th July
     8by Ole.
     9*/
     10
    511#include "Python.h"
    612#include "Numeric/arrayobject.h"
     
    2329char isdata(float);
    2430void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
    25 PyArrayObject *_read_mux2(int, char **, float *, double *, PyObject *, int);
     31float** _read_mux2(int, char **, float *, double *, int);
    2632
    2733PyObject *read_mux2(PyObject *self, PyObject *args){
     
    2935   
    3036    Python call:
    31     read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
     37    read_mux2(numSrc,filenames,weights,file_params,verbose)
    3238
    3339    NOTE:
     
    3743  PyObject *filenames;
    3844  PyArrayObject *pyweights,*file_params;
    39   PyArrayObject *sts_data;
    40   PyObject *permutation;
     45  PyArrayObject *pydata;
    4146  PyObject *fname;
    4247
     
    4651  long verbose;
    4752
     53  float **cdata;
     54  int dimensions[2];
    4855  int nsta0;
    4956  int nt;
    5057  double dt;
    51   int i;
     58  int i,j;
     59  int start_tstep;
     60  int finish_tstep;
     61  int it,time;
     62  int num_ts;
    5263 
    5364  // Convert Python arguments to C
    54   if (!PyArg_ParseTuple(args, "iOOOOi",
    55                         &numSrc,&filenames,&pyweights,&file_params,&permutation,&verbose)) {                   
     65  if (!PyArg_ParseTuple(args, "iOOOi",
     66                        &numSrc,&filenames,&pyweights,&file_params,&verbose)) {                 
    5667                       
    5768    PyErr_SetString(PyExc_RuntimeError,
     
    8394  muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *));
    8495  if (muxFileNameArray == NULL) {
    85      PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
    86      return NULL;
     96     printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
     97     exit(-1);
    8798  }
    8899  for (i=0;i<PyList_Size(filenames);i++){
    89100    muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char));
    90101    if (muxFileNameArray[i] == NULL) {
    91       PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
    92       return NULL;
     102      printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
     103      exit(-1);
    93104    }
    94105    fname=PyList_GetItem(filenames, i);
    95106    if (!PyString_Check(fname)) {
    96107      PyErr_SetString(PyExc_ValueError, "filename not a string");
    97       return NULL;
    98108    }
    99109    muxFileNameArray[i]=PyString_AsString(fname);
     
    113123
    114124  // Read in mux2 data from file
    115   sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,permutation,(int)verbose);
     125  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
     126
    116127
    117128  // Allocate space for return vector
     
    120131  nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]);
    121132
     133 
     134  // Find min and max start times of all gauges
     135  start_tstep=nt+1;
     136  finish_tstep=-1;
     137  for (i=0;i<nsta0;i++){
     138    if ((int)cdata[i][nt+3] < start_tstep){
     139      start_tstep=(int)cdata[i][nt+3];
     140    }
     141    if ((int)cdata[i][nt+4] > finish_tstep){
     142      finish_tstep=(int)cdata[i][nt+4];
     143    }
     144  }
     145
     146  if ((start_tstep>nt) | (finish_tstep < 0)){
     147    printf("ERROR: Gauge data has incorrect start and finsh times\n");
     148    return NULL;
     149  }
     150
     151  if (start_tstep>=finish_tstep){
     152    printf("ERROR: Gauge data has non-postive_length\n");
     153    return NULL;
     154  }
     155
     156  num_ts=finish_tstep-start_tstep+1;
     157  dimensions[0]=nsta0;
     158  dimensions[1]=num_ts+POFFSET;
     159  pydata = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
     160  if(pydata == NULL){
     161    printf("ERROR: Memory for pydata array could not be allocated.\n");
     162    return NULL;
     163  }
     164
     165  // Each gauge begins and ends recording at different times. When a gauge is
     166  // not recording but at least one other gauge is. Pad the non-recording gauge
     167  // array with zeros.
     168  for (i=0;i<nsta0;i++){
     169    time=0;
     170    for (it=0;it<finish_tstep;it++){
     171      if((it+1>=start_tstep)&&(it+1<=finish_tstep)){
     172        if (it+1>(int)cdata[i][nt+4]){
     173          //This gauge has stopped recording but others are still recording
     174        *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=0.0;
     175        }else{
     176          *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=cdata[i][it];
     177        }
     178        time++;
     179      }
     180    }
     181    // Pass back lat,lon,elevation
     182    for (j=0; j<POFFSET; j++){
     183      *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j];
     184     }
     185  }
     186
    122187  free(weights);
    123188  free(muxFileNameArray);
    124   return  PyArray_Return(sts_data);
     189  free(cdata);
     190  return  PyArray_Return(pydata);
    125191
    126192}
    127193
    128 PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, PyObject *permutation, int verbose)
     194float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
    129195{
    130196   FILE *fp;
    131    int nsta, nsta0, i, j, t, isrc, ista, N;
     197   int nsta, nsta0, i, isrc, ista;
    132198   struct tgsrwg *mytgs, *mytgs0;
    133199   char *muxFileName;                                                                 
     
    138204   long numData;
    139205
    140    int start_tstep;
    141    int finish_tstep;
    142    int num_ts;
    143    int dimensions[2];
     206
    144207   int len_sts_data;
     208   float **sts_data;
    145209   float *temp_sts_data;
    146 
    147    PyArrayObject *sts_data;
    148    
    149    int permuation_dimensions[1];
    150    PyArrayObject *permutation_array;
    151    
    152    long int offset;   
     210   
     211   long int offset;
     212
     213   /* Allocate space for the names and the weights and pointers to the data*/   
    153214   
    154215   /* Check that the input files have mux2 extension*/
     
    177238   /* Open the first muxfile */
    178239   if((fp=fopen(muxFileNameArray[0],"r"))==NULL){
    179      fprintf(stderr,"Cannot open file %s\n", muxFileNameArray[0]);
    180      PyErr_SetString(PyExc_RuntimeError,"");
    181      return NULL; 
     240      fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]);
     241      exit(-1); 
    182242   }
    183243 
     
    192252   /*make an array to hold the start and stop steps for each station for each
    193253     source*/   
     254   //FIXME: Should only store start and stop times for one source
    194255   fros = (int *) malloc(nsta0*numSrc*sizeof(int));
    195256   lros = (int *) malloc(nsta0*numSrc*sizeof(int));
     
    201262   /* compute the size of the data block for source 0 */   
    202263   numData = getNumData(fros, lros, nsta0);
    203    num_ts = mytgs0[0].nt;
    204    
     264
    205265   /* Burbidge: Added a sanity check here */
    206266   if (numData < 0) {
    207267     fprintf(stderr,"Size of data block appears to be negative!\n");
    208      PyErr_SetString(PyExc_RuntimeError,"");
    209      return NULL; 
     268     exit(-1);
    210269   }
    211270   fclose(fp);
    212271
    213272   // Allocate header array for each remaining source file.
     273   // FIXME: only need to store for one source and which is compared to all subsequent
     274   // source headers
    214275   if(numSrc > 1){
    215276      /* allocate space for tgsrwg for the other sources */
     
    217278   } else {
    218279     /* FIXME (JJ): What should happen in case the are no source files?*/
    219      /* If we exit here, tests will break */       
    220    }
    221    
    222    /* loop over remaining sources and read headers */
     280     /* If we exit here, tests will break */
     281     // fprintf(stderr, "No source file specified\n");
     282     // exit(-1);       
     283   }
     284   
     285   /* loop over remaining sources, check compatibility, and read them into
     286    *muxData */
    223287   for(isrc=1; isrc<numSrc; isrc++){
    224288     muxFileName = muxFileNameArray[isrc];
     
    228292       {
    229293         fprintf(stderr, "cannot open file %s\n", muxFileName);
    230          PyErr_SetString(PyExc_RuntimeError,"");
    231          return NULL;
     294         exit(-1); 
    232295       }
    233296     
     
    236299     if(nsta != nsta0){
    237300       fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
    238        PyErr_SetString(PyExc_RuntimeError,"");
    239301       fclose(fp);
    240        return NULL;
     302       exit(-1); 
    241303     }
    242304     fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
     
    245307         fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
    246308         fclose(fp);
    247          return NULL;               
     309         exit(-1);                
    248310       }   
    249        if(mytgs[ista].nt != num_ts){
     311       if(mytgs[ista].nt != mytgs0[ista].nt){
    250312         fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
    251          PyErr_SetString(PyExc_RuntimeError,"");
    252313         fclose(fp);
    253          return NULL;             
     314         exit(-1);                
    254315       }
    255316     }
     
    262323      numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0);
    263324
     325      /* Burbidge: Added a sanity check here */
    264326      if (numData < 0){
    265327          fprintf(stderr,"Size of data block appears to be negative!\n");
    266           PyErr_SetString(PyExc_RuntimeError,"");
    267           return NULL;
     328          exit(-1);
    268329      }
    269330      fclose(fp);             
    270331   }
    271 
    272    /*
    273    for (i=0;i<nsta0*numSrc;i++){
    274        printf("%d, fros %d\n",i,fros[i]);
    275        printf("%d, lros %d\n",i,lros[i]);
    276        }*/
    277 
    278    if (permutation == Py_None){
    279        // if no permutation is specified return the times series of all the gauges in the mux file
    280        permuation_dimensions[0]=nsta0;
    281        permutation_array=(PyArrayObject *) PyArray_FromDims(1, permuation_dimensions, PyArray_INT);
    282        for (ista=0; ista<nsta0; ista++){
    283            *(long *) (permutation_array -> data+ista*permutation_array->strides[0])=ista;
    284        }
    285    }else{
    286        // Specifies the gauge numbers that for which data is to be extracted
    287        permutation_array=(PyArrayObject *) PyArray_ContiguousFromObject(permutation,PyArray_INT,1,1);
    288        N = permutation_array->dimensions[0]; //FIXME: this is overwritten below
    289        for (ista=0; ista<N; ista++){
    290            if ((int)*(long *) (permutation_array -> data+ista*permutation_array->strides[0])>=nsta0){
    291                printf("Maximum index = %d, you had %d\n",nsta0-1,(int)*(long *) (permutation_array -> data+ista*permutation_array->strides[0]));
    292                PyErr_SetString(PyExc_RuntimeError,"The permutation specified is out of bounds");
    293                return NULL;
    294            }
    295        }
    296    }
    297    if(permutation_array == NULL){
    298      PyErr_SetString(PyExc_RuntimeError,"Memory for permutation_array array could not be allocated");
    299      return NULL;
    300    }
    301 
    302    N = permutation_array->dimensions[0];
    303    
    304    params[0]=(double)N;
     332   params[0]=(double)nsta0;
    305333   params[1]=(double)mytgs0[0].dt;
    306    params[2]=(double)num_ts;
    307  
    308    
    309    // Find min and max start times of all gauges
    310    start_tstep=num_ts+1;
    311    finish_tstep=-1;
    312    for (i=0;i<N;i++){
    313        ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]); 
    314        if (fros[ista]< start_tstep){
    315            start_tstep=fros[ista];
    316        }
    317        if (lros[0+nsta0*ista] > finish_tstep){
    318            finish_tstep=lros[ista];
    319        }
    320    }
    321    
    322    if ((start_tstep>num_ts) | (finish_tstep < 0)){
    323        printf("s=%d,f=%d,num_ts=%d\n",start_tstep,finish_tstep,num_ts);
    324      PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times");
    325      return NULL;
    326    }
    327    
    328    if (start_tstep>=finish_tstep){
    329      PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length");
    330      return NULL;
    331    }
    332    
    333    /* Make array(s) to hold the demuxed data */
    334    len_sts_data=num_ts+POFFSET;
    335    dimensions[0]=N;
    336    dimensions[1]=len_sts_data;
    337    sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
    338    if(sts_data == NULL){
    339     PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated");
    340     return NULL;
    341    }
    342    
    343    /* Initialise sts data to zero */
    344    for(i=0; i<N; i++){
    345      ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);
    346      for (t=0;t<num_ts;t++){
    347        *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
    348      }
    349      *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
    350      *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
    351      *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
    352      *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
    353      *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
    354    }
     334   params[2]=(double)mytgs0[0].nt;
     335   
     336   /* make array(s) to hold the demuxed data */
     337   //FIXME: This can be reduced to only contain stations in order file
     338   sts_data = (float **)malloc (nsta0 * sizeof(float *));
     339   
     340   if (sts_data == NULL){
     341     printf("ERROR: Memory for sts_data could not be allocated.\n");
     342     exit(-1);
     343   }
     344   
     345   len_sts_data=mytgs0[0].nt + POFFSET;
     346   for (ista=0; ista<nsta0; ista++){
     347     // Initialise sts_data to zero
     348     sts_data[ista] = (float *)calloc(len_sts_data, sizeof(float) );
     349     if (sts_data[ista] == NULL){
     350       printf("ERROR: Memory for sts_data could not be allocated.\n");
     351       exit(-1);
     352     }
     353     sts_data[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat;
     354     sts_data[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon;
     355     sts_data[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z;
     356     sts_data[ista][mytgs0[0].nt+3]=(float)fros[ista];
     357     sts_data[ista][mytgs0[0].nt+4]=(float)lros[ista];
     358   }
    355359
    356360   temp_sts_data = (float *)calloc(len_sts_data, sizeof(float) );
     
    365369     if((fp=fopen(muxFileName,"r"))==NULL)
    366370       {
    367          fprintf(stderr, "Cannot open file %s\n", muxFileName);
    368          PyErr_SetString(PyExc_RuntimeError,"");
    369          return NULL;
     371         //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
     372         fprintf(stderr, "cannot open file %s\n", muxFileName);
     373         exit(-1); 
    370374       }
    371375       
     
    383387
    384388     /* loop over stations */
    385      for(j=0; j<N; j++){
    386        ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]);
    387        //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0]));
     389     for(ista=0; ista<nsta0; ista++){               
    388390       
    389391       /* fill the data0 array from the mux file, and weight it */
    390        fillDataArray(ista, nsta0, num_ts, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
     392       fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
    391393       
    392394       /* weight appropriately and add */
    393        for(t=0; t<num_ts; t++){
    394          if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){
    395            *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc];
     395       for(i=0; i<mytgs0[ista].nt; i++){
     396         if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])){
     397           sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
     398           //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
    396399         }else{
    397            *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0;
     400           sts_data[ista][i] = NODATA;
    398401         }
    399402       }
    400403     }
    401    }   
     404   }
    402405   free(muxData);
    403406   free(temp_sts_data);
    404    free(fros);
    405    free(lros);
    406407   return sts_data;
    407408}   
     
    484485}
    485486
    486  
     487/* Burbidge: No "offset" is sent. Replace with max. Added grid_id */
     488void wrttxt( fname, dt, nt, x, beg, lat, lon, max, depth, grid_id, istart, istop )
     489char *fname;
     490float dt, *x, beg, max, lat, lon, depth;
     491int grid_id;
     492int nt;
     493int istart, istop;
     494{
     495   int it;
     496   float t;
     497   FILE *fp;
     498
     499   fp = fopen(fname,"w");
     500   fprintf(fp,"> lat %g lon %g depth %g max %g start_time %g stop_time %g grid_id %d\n", lat, lon, depth, max,
     501            istart*dt, istop*dt, grid_id );
     502   for(it=0;it<nt;it++)
     503   {
     504      t=beg+it*dt;
     505      fprintf(fp,"%9.3f %g\n", t, x[it]);
     506   }
     507   fclose(fp);
     508}
     509 
    487510char isdata(float x)
    488511{
Note: See TracChangeset for help on using the changeset viewer.