Changeset 5485 for anuga_core/source

Timestamp:

Jul 10, 2008, 2:14:22 PM (17 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:

• data_manager.py (modified) (7 diffs)
• test_data_manager.py (modified) (5 diffs)
• urs_ext.c (modified) (20 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -48 +48 @@
 
 """
+
+# This file was reverted from changeset:5484 to changeset:5470 on 10th July 
+# by Ole.
 
 import exceptions
@@ -4838 +4841 @@
 NORTH_VELOCITY_MUX2_LABEL =  '-n-mux2'    
 
-def read_mux2_py(filenames,weights,permutation=None,verbose=False):
-    """
-    Access the mux files specified in the filenames list. Combine the
-    data found therin as a weighted linear sum as specifed by the weights.
-    If permutation is None extract timeseries data for all gauges within the
-    files.
-    
-    Input:
-    filenames:   List of filenames specifiying the file containing the
-                 timeseries data (mux2 format) for each source
-    weights:     Weights associated with each source
-    permutation: Specifies the gauge numbers that for which data is to be
-                 extracted
-    """
+def read_mux2_py(filenames,weights,verbose=False):
+
     from Numeric import ones,Float,compress,zeros,arange
     from urs_ext import read_mux2
@@ -4865 +4856 @@
         verbose=0
         
-    data=read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
+    data=read_mux2(numSrc,filenames,weights,file_params,verbose)
 
     msg='File parameter values were not read in correctly from c file'
@@ -5046 +5037 @@
                 raise IOError, msg
 
-    if ordering_filename is not None:
-        # Read ordering file 
-        
-        try:
-            fid=open(ordering_filename, 'r')
-            file_header=fid.readline().split(',')
-            ordering_lines=fid.readlines()
-            fid.close()
-        except:
-            msg = 'Cannot open %s'%ordering_filename
-            raise Exception, msg
-
-        reference_header = 'index, longitude, latitude\n'
-        reference_header_split = reference_header.split(',')
-        for i in range(3):
-            if not file_header[i].strip() == reference_header_split[i].strip():
-                msg = 'File must contain header: '+reference_header+'\n'
-                raise Exception, msg
-
-        if len(ordering_lines)<2:
-            msg = 'File must contain at least two points'
-            raise Exception, msg
-
-   
-        
-        permutation = [int(line.split(',')[0]) for line in ordering_lines]
-    else:
-        permutation = None  
-
+    #need to do this for velocity-e-mux2 and velocity-n-mux2 files as well
+    #for now set x_momentum and y_moentum quantities to zero
     if (verbose): print 'reading mux2 file'
     mux={}
     for i,quantity in enumerate(quantities):
         # For each quantity read the associated list of source mux2 file with extenstion associated with that quantity
-        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,permutation=permutation,verbose=verbose)
+        times, latitudes_urs, longitudes_urs, elevation, mux[quantity],starttime = read_mux2_py(files_in[i],weights,verbose=verbose)
         if quantity!=quantities[0]:
             msg='%s, %s and %s have inconsitent gauge data'%(files_in[0],files_in[1],files_in[2])
@@ -5092 +5056 @@
         elevation_old=elevation
         starttime_old=starttime
+
         
     if ordering_filename is not None:
-        latitudes = latitudes_urs
-        longitudes = longitudes_urs
+        # Read ordering file 
+        
+        try:
+            fid=open(ordering_filename, 'r')
+            file_header=fid.readline().split(',')
+            ordering_lines=fid.readlines()
+            fid.close()
+        except:
+            msg = 'Cannot open %s'%ordering_filename
+            raise Exception, msg
+
+        reference_header = 'index, longitude, latitude\n'
+        reference_header_split = reference_header.split(',')
+        for i in range(3):
+            if not file_header[i] == reference_header_split[i]:
+                msg = 'File must contain header: '+reference_header+'\n'
+                raise Exception, msg
+
+        if len(ordering_lines)<2:
+            msg = 'File must contain at least two points'
+            raise Exception, msg
+
+   
+        
+        permutation = [int(line.split(',')[0]) for line in ordering_lines]
+   
+        latitudes = take(latitudes_urs, permutation)
+        longitudes = take(longitudes_urs, permutation)
         
         # Self check - can be removed to improve speed
@@ -5128 +5119 @@
         longitudes=longitudes_urs
         permutation = range(latitudes.shape[0])                
+        
 
     msg='File is empty and or clipped region not in file region'
@@ -5189 +5181 @@
     if verbose: print 'Converting quantities'
     for j in range(len(times)):
-        for i in range(number_of_points):
-            w = zscale*mux['HA'][i,j] + mean_stage
-            h=w-elevation[i]
+        for i, index in enumerate(permutation):
+            w = zscale*mux['HA'][index,j] + mean_stage
+            h=w-elevation[index]
             stage[j,i] = w
 
-            xmomentum[j,i] = mux['UA'][i,j]*h
-            ymomentum[j,i] = mux['VA'][i,j]*h
+            xmomentum[j,i] = mux['UA'][index,j]*h
+            ymomentum[j,i] = mux['VA'][index,j]*h
 
     outfile.close()

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -2 +2 @@
 #
 
+# This file was reverted from changeset:5484 to changeset:5470 on 10th July 
+# by Ole.
 
 import unittest
@@ -6105 +6107 @@
         return base_name, files
 
-    def test_read_mux2_py0(self):
+    def test_read_mux2_pyI(self):
         """Constant stage,momentum at each gauge
         """
@@ -6543 +6545 @@
         """Test multiple sources with ordering file
         """
+        
         tide=0
         time_step_count = 3
@@ -6944 +6947 @@
             fid.write(line)
         fid.close()
-        
+
         sts_file=base_name
         urs2sts(base_name, basename_out=sts_file,
@@ -8888 +8891 @@
 if __name__ == "__main__":
 
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_read_mux2_py0')
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -3 +3 @@
 gcc -shared urs_ext.o  -o urs_ext.so
 */
+
+/*
+This file was reverted from changeset:5484 to changeset:5470 on 10th July 
+by Ole.
+*/
+
 #include "Python.h"
 #include "Numeric/arrayobject.h"
@@ -23 +29 @@
 char isdata(float);
 void wrttxt(char *, float, int, float *, float, float, float, float, float, int, int, int);
-PyArrayObject *_read_mux2(int, char **, float *, double *, PyObject *, int);
+float** _read_mux2(int, char **, float *, double *, int);
 
 PyObject *read_mux2(PyObject *self, PyObject *args){
@@ -29 +35 @@
    
     Python call:
-    read_mux2(numSrc,filenames,weights,file_params,permutation,verbose)
+    read_mux2(numSrc,filenames,weights,file_params,verbose)
 
     NOTE:
@@ -37 +43 @@
   PyObject *filenames;
   PyArrayObject *pyweights,*file_params;
-  PyArrayObject *sts_data;
-  PyObject *permutation;
+  PyArrayObject *pydata;
   PyObject *fname;
 
@@ -46 +51 @@
   long verbose;
 
+  float **cdata;
+  int dimensions[2];
   int nsta0;
   int nt;
   double dt;
-  int i;
+  int i,j;
+  int start_tstep;
+  int finish_tstep;
+  int it,time;
+  int num_ts;
   
   // Convert Python arguments to C
-  if (!PyArg_ParseTuple(args, "iOOOOi",
-                        &numSrc,&filenames,&pyweights,&file_params,&permutation,&verbose)) {                    
+  if (!PyArg_ParseTuple(args, "iOOOi",
+                        &numSrc,&filenames,&pyweights,&file_params,&verbose)) {                 
                         
     PyErr_SetString(PyExc_RuntimeError, 
@@ -83 +94 @@
   muxFileNameArray = (char **) malloc((int)numSrc*sizeof(char *));
   if (muxFileNameArray == NULL) {
-     PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
-     return NULL;
+     printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
+     exit(-1);
   }
   for (i=0;i<PyList_Size(filenames);i++){
     muxFileNameArray[i] = (char *) malloc((MAX_FILE_NAME_LENGTH+1)*sizeof(char));
     if (muxFileNameArray[i] == NULL) {
-      PyErr_SetString(PyExc_RuntimeError,"Memory for muxFileNameArray could not be allocated");
-      return NULL;
+      printf("ERROR: Memory for muxFileNameArray could not be allocated.\n");
+      exit(-1);
     }
     fname=PyList_GetItem(filenames, i);
     if (!PyString_Check(fname)) {
       PyErr_SetString(PyExc_ValueError, "filename not a string");
-      return NULL;
     }
     muxFileNameArray[i]=PyString_AsString(fname);
@@ -113 +123 @@
 
   // Read in mux2 data from file
-  sts_data=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,permutation,(int)verbose);
+  cdata=_read_mux2((int)numSrc,muxFileNameArray,weights,(double*)file_params->data,(int)verbose);
+
 
   // Allocate space for return vector
@@ -120 +131 @@
   nt=(int)*(double *) (file_params -> data+2*file_params->strides[0]);
 
+  
+  // Find min and max start times of all gauges
+  start_tstep=nt+1;
+  finish_tstep=-1;
+  for (i=0;i<nsta0;i++){
+    if ((int)cdata[i][nt+3] < start_tstep){
+      start_tstep=(int)cdata[i][nt+3];
+    }
+    if ((int)cdata[i][nt+4] > finish_tstep){
+      finish_tstep=(int)cdata[i][nt+4]; 
+    }
+  }
+
+  if ((start_tstep>nt) | (finish_tstep < 0)){
+    printf("ERROR: Gauge data has incorrect start and finsh times\n");
+    return NULL;
+  }
+
+  if (start_tstep>=finish_tstep){
+    printf("ERROR: Gauge data has non-postive_length\n");
+    return NULL;
+  }
+
+  num_ts=finish_tstep-start_tstep+1;
+  dimensions[0]=nsta0;
+  dimensions[1]=num_ts+POFFSET;
+  pydata = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
+  if(pydata == NULL){
+    printf("ERROR: Memory for pydata array could not be allocated.\n");
+    return NULL;
+  }
+
+  // Each gauge begins and ends recording at different times. When a gauge is
+  // not recording but at least one other gauge is. Pad the non-recording gauge
+  // array with zeros.
+  for (i=0;i<nsta0;i++){
+    time=0;
+    for (it=0;it<finish_tstep;it++){
+      if((it+1>=start_tstep)&&(it+1<=finish_tstep)){
+        if (it+1>(int)cdata[i][nt+4]){
+          //This gauge has stopped recording but others are still recording
+        *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=0.0;
+        }else{
+          *(double *) (pydata -> data+i*pydata->strides[0]+time*pydata->strides[1])=cdata[i][it];
+        }
+        time++;
+      }
+    }
+    // Pass back lat,lon,elevation
+    for (j=0; j<POFFSET; j++){
+      *(double *) (pydata -> data+i*pydata->strides[0]+(num_ts+j)*pydata->strides[1])=cdata[i][nt+j];
+     }
+  }
+
   free(weights);
   free(muxFileNameArray);
-  return  PyArray_Return(sts_data);
+  free(cdata);
+  return  PyArray_Return(pydata);
 
 }
 
-PyArrayObject *_read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, PyObject *permutation, int verbose)
+float** _read_mux2(int numSrc, char **muxFileNameArray, float *weights, double *params, int verbose)
 {
    FILE *fp;
-   int nsta, nsta0, i, j, t, isrc, ista, N;
+   int nsta, nsta0, i, isrc, ista;
    struct tgsrwg *mytgs, *mytgs0;
    char *muxFileName;                                                                  
@@ -138 +204 @@
    long numData;
 
-   int start_tstep;
-   int finish_tstep;
-   int num_ts;
-   int dimensions[2];
+
    int len_sts_data;
+   float **sts_data;
    float *temp_sts_data;
-
-   PyArrayObject *sts_data;
-   
-   int permuation_dimensions[1];
-   PyArrayObject *permutation_array;
-   
-   long int offset;   
+   
+   long int offset;
+
+   /* Allocate space for the names and the weights and pointers to the data*/    
    
    /* Check that the input files have mux2 extension*/
@@ -177 +238 @@
    /* Open the first muxfile */
    if((fp=fopen(muxFileNameArray[0],"r"))==NULL){
-     fprintf(stderr,"Cannot open file %s\n", muxFileNameArray[0]);
-     PyErr_SetString(PyExc_RuntimeError,"");
-     return NULL;  
+      fprintf(stderr, "cannot open file %s\n", muxFileNameArray[0]);
+      exit(-1);  
    }
  
@@ -192 +252 @@
    /*make an array to hold the start and stop steps for each station for each
      source*/   
+   //FIXME: Should only store start and stop times for one source
    fros = (int *) malloc(nsta0*numSrc*sizeof(int));
    lros = (int *) malloc(nsta0*numSrc*sizeof(int));
@@ -201 +262 @@
    /* compute the size of the data block for source 0 */   
    numData = getNumData(fros, lros, nsta0);
-   num_ts = mytgs0[0].nt;
-   
+
    /* Burbidge: Added a sanity check here */
    if (numData < 0) {
      fprintf(stderr,"Size of data block appears to be negative!\n");
-     PyErr_SetString(PyExc_RuntimeError,"");
-     return NULL;  
+     exit(-1);
    }
    fclose(fp); 
 
    // Allocate header array for each remaining source file.
+   // FIXME: only need to store for one source and which is compared to all subsequent
+   // source headers
    if(numSrc > 1){
       /* allocate space for tgsrwg for the other sources */
@@ -217 +278 @@
    } else {
      /* FIXME (JJ): What should happen in case the are no source files?*/
-     /* If we exit here, tests will break */       
-   }
-   
-   /* loop over remaining sources and read headers */
+     /* If we exit here, tests will break */
+     // fprintf(stderr, "No source file specified\n");
+     // exit(-1);       
+   }
+   
+   /* loop over remaining sources, check compatibility, and read them into 
+    *muxData */
    for(isrc=1; isrc<numSrc; isrc++){
      muxFileName = muxFileNameArray[isrc];
@@ -228 +292 @@
        {
          fprintf(stderr, "cannot open file %s\n", muxFileName);
-         PyErr_SetString(PyExc_RuntimeError,"");
-         return NULL; 
+         exit(-1);  
        }
      
@@ -236 +299 @@
      if(nsta != nsta0){
        fprintf(stderr,"%s has different number of stations to %s\n", muxFileName, muxFileNameArray[0]);
-       PyErr_SetString(PyExc_RuntimeError,"");
        fclose(fp);
-       return NULL; 
+       exit(-1);   
      }
      fread(&mytgs[0], nsta*sizeof(struct tgsrwg), 1, fp);
@@ -245 +307 @@
          fprintf(stderr,"%s has different sampling rate to %s\n", muxFileName, muxFileNameArray[0]);
          fclose(fp);
-         return NULL;                
+         exit(-1);                 
        }   
-       if(mytgs[ista].nt != num_ts){
+       if(mytgs[ista].nt != mytgs0[ista].nt){
          fprintf(stderr,"%s has different series length to %s\n", muxFileName, muxFileNameArray[0]);
-         PyErr_SetString(PyExc_RuntimeError,"");
          fclose(fp);
-         return NULL;              
+         exit(-1);                 
        }
      }
@@ -262 +323 @@
       numData = getNumData(fros+isrc*nsta0, lros+isrc*nsta0, nsta0);
 
+      /* Burbidge: Added a sanity check here */
       if (numData < 0){
           fprintf(stderr,"Size of data block appears to be negative!\n");
-          PyErr_SetString(PyExc_RuntimeError,"");
-          return NULL;
+          exit(-1);
       }
       fclose(fp);             
    }
-
-   /*
-   for (i=0;i<nsta0*numSrc;i++){
-       printf("%d, fros %d\n",i,fros[i]);
-       printf("%d, lros %d\n",i,lros[i]);
-       }*/
-
-   if (permutation == Py_None){
-       // if no permutation is specified return the times series of all the gauges in the mux file
-       permuation_dimensions[0]=nsta0;
-       permutation_array=(PyArrayObject *) PyArray_FromDims(1, permuation_dimensions, PyArray_INT);
-       for (ista=0; ista<nsta0; ista++){
-           *(long *) (permutation_array -> data+ista*permutation_array->strides[0])=ista;
-       } 
-   }else{
-       // Specifies the gauge numbers that for which data is to be extracted
-       permutation_array=(PyArrayObject *) PyArray_ContiguousFromObject(permutation,PyArray_INT,1,1);
-       N = permutation_array->dimensions[0]; //FIXME: this is overwritten below
-       for (ista=0; ista<N; ista++){
-           if ((int)*(long *) (permutation_array -> data+ista*permutation_array->strides[0])>=nsta0){
-               printf("Maximum index = %d, you had %d\n",nsta0-1,(int)*(long *) (permutation_array -> data+ista*permutation_array->strides[0]));
-               PyErr_SetString(PyExc_RuntimeError,"The permutation specified is out of bounds");
-               return NULL;
-           }
-       }
-   }
-   if(permutation_array == NULL){
-     PyErr_SetString(PyExc_RuntimeError,"Memory for permutation_array array could not be allocated");
-     return NULL;
-   }
-
-   N = permutation_array->dimensions[0];
-   
-   params[0]=(double)N;
+   params[0]=(double)nsta0;
    params[1]=(double)mytgs0[0].dt;
-   params[2]=(double)num_ts;
-  
-    
-   // Find min and max start times of all gauges
-   start_tstep=num_ts+1;
-   finish_tstep=-1;
-   for (i=0;i<N;i++){
-       ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);  
-       if (fros[ista]< start_tstep){
-           start_tstep=fros[ista];
-       }
-       if (lros[0+nsta0*ista] > finish_tstep){
-           finish_tstep=lros[ista]; 
-       }
-   }
-   
-   if ((start_tstep>num_ts) | (finish_tstep < 0)){
-       printf("s=%d,f=%d,num_ts=%d\n",start_tstep,finish_tstep,num_ts);
-     PyErr_SetString(PyExc_RuntimeError,"Gauge data has incorrect start and finsh times");
-     return NULL;
-   }
-   
-   if (start_tstep>=finish_tstep){
-     PyErr_SetString(PyExc_RuntimeError,"Gauge data has non-postive_length");
-     return NULL;
-   }
-   
-   /* Make array(s) to hold the demuxed data */
-   len_sts_data=num_ts+POFFSET;
-   dimensions[0]=N;
-   dimensions[1]=len_sts_data;
-   sts_data = (PyArrayObject *) PyArray_FromDims(2, dimensions, PyArray_DOUBLE);
-   if(sts_data == NULL){
-    PyErr_SetString(PyExc_RuntimeError,"Memory for pydata array could not be allocated");
-    return NULL;
-   }
-   
-   /* Initialise sts data to zero */
-   for(i=0; i<N; i++){ 
-     ista = *(long *) (permutation_array -> data+i*permutation_array->strides[0]);
-     for (t=0;t<num_ts;t++){
-       *(double *) (sts_data -> data+i*sts_data->strides[0]+t*sts_data->strides[1])=0.0;
-     }
-     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts)*sts_data->strides[1])=(float)mytgs0[ista].geolat;
-     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+1)*sts_data->strides[1])=(float)mytgs0[ista].geolon;
-     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+2)*sts_data->strides[1])=(float)mytgs0[ista].z;
-     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+3)*sts_data->strides[1])=(float)fros[ista];
-     *(double *) (sts_data -> data+i*sts_data->strides[0]+(num_ts+4)*sts_data->strides[1])=(float)lros[ista];
-   } 
+   params[2]=(double)mytgs0[0].nt;
+   
+   /* make array(s) to hold the demuxed data */
+   //FIXME: This can be reduced to only contain stations in order file
+   sts_data = (float **)malloc (nsta0 * sizeof(float *));
+   
+   if (sts_data == NULL){
+     printf("ERROR: Memory for sts_data could not be allocated.\n");
+     exit(-1);
+   }
+   
+   len_sts_data=mytgs0[0].nt + POFFSET;
+   for (ista=0; ista<nsta0; ista++){
+     // Initialise sts_data to zero
+     sts_data[ista] = (float *)calloc(len_sts_data, sizeof(float) );
+     if (sts_data[ista] == NULL){
+       printf("ERROR: Memory for sts_data could not be allocated.\n");
+       exit(-1);
+     }
+     sts_data[ista][mytgs0[0].nt]=(float)mytgs0[ista].geolat;
+     sts_data[ista][mytgs0[0].nt+1]=(float)mytgs0[ista].geolon;
+     sts_data[ista][mytgs0[0].nt+2]=(float)mytgs0[ista].z;
+     sts_data[ista][mytgs0[0].nt+3]=(float)fros[ista];
+     sts_data[ista][mytgs0[0].nt+4]=(float)lros[ista];
+   }
 
    temp_sts_data = (float *)calloc(len_sts_data, sizeof(float) );
@@ -365 +369 @@
      if((fp=fopen(muxFileName,"r"))==NULL)
        {
-         fprintf(stderr, "Cannot open file %s\n", muxFileName);
-         PyErr_SetString(PyExc_RuntimeError,"");
-         return NULL;
+         //fprintf(stderr, "%s: cannot open file %s\n", av[0], muxFileName);
+         fprintf(stderr, "cannot open file %s\n", muxFileName);
+         exit(-1);  
        }
        
@@ -383 +387 @@
 
      /* loop over stations */
-     for(j=0; j<N; j++){ 
-       ista = *(long *) (permutation_array -> data+j*permutation_array->strides[0]);
-       //printf("%d\n",(int)*(long *) (permutation_array -> data+j*permutation_array->strides[0]));
+     for(ista=0; ista<nsta0; ista++){               
        
        /* fill the data0 array from the mux file, and weight it */
-       fillDataArray(ista, nsta0, num_ts, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
+       fillDataArray(ista, nsta0, mytgs0[ista].nt, mytgs0[ista].ig, fros+isrc*nsta0, lros+isrc*nsta0, temp_sts_data, &istart, &istop, muxData);
        
        /* weight appropriately and add */
-       for(t=0; t<num_ts; t++){
-         if((isdata(*(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]))) && isdata(temp_sts_data[t])){
-           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) += temp_sts_data[t] * weights[isrc];
+       for(i=0; i<mytgs0[ista].nt; i++){
+         if((isdata(sts_data[ista][i])) && isdata(temp_sts_data[i])){
+           sts_data[ista][i] += temp_sts_data[i] * weights[isrc];
+           //printf("%d,%d,%f\n",ista,i,sts_data[ista][i]);
          }else{
-           *(double *) (sts_data -> data+j*sts_data->strides[0]+t*sts_data->strides[1]) = 0.0;
+           sts_data[ista][i] = NODATA;
          }
        }
      }
-   }   
+   }
    free(muxData);
    free(temp_sts_data);
-   free(fros);
-   free(lros);
    return sts_data;
 }   
@@ -484 +485 @@
 } 
 
- 
+/* Burbidge: No "offset" is sent. Replace with max. Added grid_id */
+void wrttxt( fname, dt, nt, x, beg, lat, lon, max, depth, grid_id, istart, istop )
+char *fname;
+float dt, *x, beg, max, lat, lon, depth;
+int grid_id;
+int nt;
+int istart, istop;
+{
+   int it;
+   float t;
+   FILE *fp;
+
+   fp = fopen(fname,"w");
+   fprintf(fp,"> lat %g lon %g depth %g max %g start_time %g stop_time %g grid_id %d\n", lat, lon, depth, max, 
+            istart*dt, istop*dt, grid_id );
+   for(it=0;it<nt;it++)
+   {
+      t=beg+it*dt;
+      fprintf(fp,"%9.3f %g\n", t, x[it]);
+   }
+   fclose(fp);
+}
+  
 char isdata(float x)
 {