Changeset 5560 for anuga_core/source

Timestamp:

Jul 23, 2008, 1:38:52 PM (17 years ago)

Author:

nariman

Message:

Creation of static work arrays. Cleanup of code.

File:

• anuga_core/source/anuga/shallow_water/urs_ext.c (modified) (4 diffs)

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -25 +25 @@
 #define POFFSET 5 //Number of site_params
 
-void fillDataArray(int, int, int, int, int *, int *, float *, 
-           int *, int *, float *);
-long getNumData(const int *fros, const int *lros, const int nsta);
-char isdata(float);
-void _read_mux2_headers(int numSrc, 
-                           char **muxFileNameArray, 
-                   double *params,
-                                   int** fros, int** lros, struct tgsrwg ** mytgs0,
-                   int* nsta0,
-                   int verbose);
+static int *fros=NULL; 
+static int *lros=NULL;
+static struct tgsrwg* mytgs0=NULL;
+
+static long numDataMax=0;
+
+/////////////////////////////////////////////////////////////////////////
+//Auxiliary functions
+void fillDataArray(int ista, int total_number_of_stations, int nt, int ig, int *nst, 
+                   int *nft, float *data, int *istart_p, 
+           int *istop_p, float *muxData)
+{
+    int it, last_it, jsta;
+    long int offset=0;
+
+
+    last_it = -1;
+    /* make arrays of starting and finishing time steps for the tide gauges */
+    /* and fill them from the file */
+
+    /* update start and stop timesteps for this gauge */
+    if (nst[ista]!= -1)
+    {
+        if(*istart_p == -1)
+        {
+            *istart_p = nst[ista];
+        }
+        else
+        {
+            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
+        }
+    }
+    if (nft[ista] != -1)
+    {
+        if (*istop_p == -1)
+        {
+            *istop_p = nft[ista];
+        }
+        else
+        {
+            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
+        }
+    }     
+    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
+    {
+        /* gauge never started recording, or was outside of all grids, 
+    fill array with 0 */
+        for(it = 0; it < nt; it++)
+        {
+            data[it] = 0.0;
+        }
+    }   
+    else
+    {
+        for(it = 0; it < nt; it++)
+        {
+            last_it = it;
+            /* skip t record of data block */
+            offset++;
+            /* skip records from earlier tide gauges */
+            for(jsta = 0; jsta < ista; jsta++)
+                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
+                    offset++;
+
+            /* deal with the tide gauge at hand */
+            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
+                /* gauge is recording at this time */
+            {
+                memcpy(data + it, muxData + offset, sizeof(float));
+                offset++;
+            }
+            else if (it + 1 < nst[ista])
+            {
+                /* gauge has not yet started recording */
+                data[it] = 0.0;
+            }   
+            else
+                /* gauge has finished recording */
+            {
+                data[it] = NODATA;
+                break;
+            }
+
+            /* skip records from later tide gauges */
+            for(jsta = ista + 1; jsta < total_number_of_stations; jsta++)
+                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
+                    offset++;
+        }
+
+        if(last_it < nt - 1)
+            /* the loop was exited early because the gauge had 
+        finished recording */
+            for(it = last_it+1; it < nt; it++)
+                data[it] = NODATA;
+    }
+} 
+
+
+char isdata(float x)
+{
+    //char value;
+    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
+    {
+       return 0;
+    }
+    else
+    {
+        return 1;  
+    }
+}
+
+
+long getNumData(const int *fros, const int *lros, const int total_number_of_stations)
+/* calculates the number of data in the data block of a mux file */
+/* based on the first and last recorded output steps for each gauge */ 
+{
+    int ista, last_output_step;
+    long numData = 0;
+
+    last_output_step = 0;   
+    for(ista = 0; ista < total_number_of_stations; ista++)
+        if(*(fros + ista) != -1)
+        {
+            numData += *(lros + ista) - *(fros + ista) + 1;
+            last_output_step = (last_output_step < *(lros+ista) ? 
+                            *(lros+ista):last_output_step);
+        }   
+        numData += last_output_step*total_number_of_stations; /* these are the t records */
+        return numData;
+}
+
+/////////////////////////////////////////////////////////////////////////
+//Internal Functions
+int _read_mux2_headers(int numSrc, 
+                        char **muxFileNameArray, 
+                        int* total_number_of_stations,
+                                                int* number_of_time_steps,
+                                                double* delta_t,
+                        //long* numDataMax,
+                        int verbose)
+{
+    FILE *fp;
+    int numsta, i, j;
+    struct tgsrwg *mytgs=0;
+    char *muxFileName;                                                                  
+    char susMuxFileName;
+    long numData;
+
+    /* Allocate space for the names and the weights and pointers to the data*/
+
+    /* Check that the input files have mux2 extension*/
+    susMuxFileName = 0;
+    for(i = 0; i < numSrc; i++)
+    { 
+        muxFileName = muxFileNameArray[i];
+        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 
+                     "mux2") != 0)
+        {
+            susMuxFileName = 1;
+            break;
+        }
+    }
+
+    if(susMuxFileName)
+    {
+        printf("\n**************************************************************************\n");
+        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
+        printf("   At least one input file name does not end with mux2\n");
+        printf("   Check your results carefully!\n");
+        printf("**************************************************************************\n\n");
+    }   
+
+    if (verbose)
+    {
+        printf("Reading mux header information\n");
+    }
+
+    /* Loop over all sources, read headers and check compatibility */
+    for (i = 0; i < numSrc; i++)
+    {
+        muxFileName = muxFileNameArray[i];
+
+        /* open the mux file */
+        if((fp = fopen(muxFileName, "r")) == NULL)
+        {
+            fprintf(stderr, "cannot open file %s\n", muxFileName);
+            return 0;  
+        }
+        
+        if (!i)
+        {
+            fread(total_number_of_stations, sizeof(int), 1, fp);
+        
+            fros = (int*)malloc(*total_number_of_stations*numSrc*sizeof(int)); 
+            lros = (int*)malloc(*total_number_of_stations*numSrc*sizeof(int));
+      
+            mytgs0 = (struct tgsrwg*)malloc(*total_number_of_stations*sizeof(struct tgsrwg));
+            mytgs = (struct tgsrwg*)malloc(*total_number_of_stations*sizeof(struct tgsrwg));
+
+            fread(mytgs0, *total_number_of_stations*sizeof(struct tgsrwg), 1, fp);
+        }
+        else
+        {
+            /* check that the mux files are compatible */      
+            fread(&numsta, sizeof(int), 1, fp);
+            if(numsta != *total_number_of_stations)
+            {
+                fprintf(stderr,"%s has different number of stations to %s\n", 
+                muxFileName, 
+                muxFileNameArray[0]);
+                fclose(fp);
+                return 0;   
+            }
+
+            fread(mytgs, numsta*sizeof(struct tgsrwg), 1, fp); 
+            
+            for (j = 0; j < numsta; j++)
+            {
+                if (mytgs[j].dt != mytgs0[j].dt)
+                {
+                    fprintf(stderr,"%s has different sampling rate to %s\n", 
+                    muxFileName, 
+                    muxFileNameArray[0]);
+                    fclose(fp);
+                    return 0;            
+                }   
+                if (mytgs[j].nt != mytgs0[j].nt)
+                {
+                    fprintf(stderr,"%s has different series length to %s\n", 
+                    muxFileName, 
+                    muxFileNameArray[0]);
+                    fclose(fp);
+                    return 0;            
+                }
+
+                if (mytgs[j].nt != mytgs0[0].nt)
+                {
+                    printf("Station 0 has different series length to Station %d\n", j); 
+                }
+            }
+        }
+
+        /* Read the start and stop times for this source */
+        fread(fros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
+        fread(lros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
+
+        /* Compute the size of the data block for this source */
+        numData = getNumData(fros + i*(*total_number_of_stations), lros + i*(*total_number_of_stations), (*total_number_of_stations));
+
+        /* Sanity check */
+        if (numData < 0)
+        {
+            fprintf(stderr,"Size of data block appears to be negative!\n");
+            return 0;        
+        }
+
+        if (numDataMax < numData)
+        {
+            numDataMax = numData;
+        }
+
+        fclose(fp);          
+    }
+
+    *delta_t = (double)mytgs0[0].dt;
+    *number_of_time_steps = mytgs0[0].nt;
+
+        free(mytgs);
+
+    return 1;
+}
+
 
 float** _read_mux2(int numSrc, 
-           char **muxFileNameArray, 
-           float *weights, 
-           double *params, 
-           int *number_of_stations,
-           long *permutation,
-           int verbose);
-
-    PyObject *read_mux2(PyObject *self, PyObject *args){
+                   char **muxFileNameArray, 
+                   float *weights, 
+                   double *params, 
+                   int *number_of_stations,
+                   long *permutation,
+                   int verbose)
+{
+    FILE *fp;
+    int total_number_of_stations, i, isrc, ista, k;
+    //struct tgsrwg* mytgs0=NULL;
+    char *muxFileName;
+    int istart, istop;
+    int number_of_selected_stations;
+    float *muxData=NULL; // Suppress warning
+    long numData;
+
+    int len_sts_data;
+    float **sts_data;
+    float *temp_sts_data;
+
+    long int offset;
+
+    int number_of_time_steps;
+    double delta_t;
+    //long numDataMax;
+    
+    _read_mux2_headers(numSrc, 
+                       muxFileNameArray, 
+                       &total_number_of_stations,
+                       &number_of_time_steps,
+                       &delta_t,
+                       verbose);
+
+    //FIXME: The params can be removed after the python interface is modified.
+    params[0] = (double)total_number_of_stations;
+    params[1] = (double)delta_t;
+    params[2] = (double)number_of_time_steps;
+    
+    // Apply rule that an empty permutation file means 'take all stations'
+    // We could change this later by passing in None instead of the empty 
+    // permutation.
+    number_of_selected_stations = *number_of_stations;  
+    if (number_of_selected_stations == 0)
+    {
+        number_of_selected_stations = total_number_of_stations;  
+    
+        // Return possibly updated number of stations
+        *number_of_stations = total_number_of_stations;     
+      
+        // Create the Identity permutation vector
+        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
+        for (i = 0; i < number_of_selected_stations; i++)
+        {
+            permutation[i] = (long) i;  
+        }
+    }
+    
+    // Make array(s) to hold demuxed data for stations given in the 
+    // permutation file 
+    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
+    if (sts_data == NULL)
+    {
+        printf("ERROR: Memory for sts_data could not be allocated.\n");
+        return NULL;
+    }
+
+    // For each selected station, allocate space for its data
+
+    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
+    for (i = 0; i < number_of_selected_stations; i++)
+    {
+        // Initialise sts_data to zero
+        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
+        if (sts_data[i] == NULL)
+        {
+            printf("ERROR: Memory for sts_data could not be allocated.\n");
+            return NULL;
+        }
+
+        ista = (int) permutation[i]; // Get global index into mux data
+    
+        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
+        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
+        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
+        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
+        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
+    }
+
+    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
+
+    muxData = (float*)malloc(numDataMax*sizeof(float));
+    /* Loop over all sources */
+    //FIXME: remove istart and istop they are not used.
+    istart = -1;
+    istop = -1;
+    for (isrc = 0; isrc < numSrc; isrc++)
+    {
+        /* Read in data block from mux2 file */
+        muxFileName = muxFileNameArray[isrc];
+        if((fp = fopen(muxFileName, "r")) == NULL)
+        {
+            fprintf(stderr, "cannot open file %s\n", muxFileName);
+            return NULL;                    
+        }
+
+        if (verbose){
+            printf("Reading mux file %s\n", muxFileName);
+        }
+
+        offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
+        fseek(fp, offset, 0);
+
+        numData = getNumData(fros + isrc*total_number_of_stations, lros + isrc*total_number_of_stations, total_number_of_stations);
+        fread(muxData, numData*sizeof(float), 1, fp); 
+        fclose(fp);
+
+        // loop over stations present in the permutation array 
+        //     use ista with mux data
+        //     use i with the processed data to be returned         
+        
+        for(i = 0; i < number_of_selected_stations; i++)
+        {               
+    
+            ista = (int) permutation[i]; // Get global index into mux data  
+        
+            /* fill the data0 array from the mux file, and weight it */
+            fillDataArray(ista, 
+                          total_number_of_stations, 
+                          mytgs0[ista].nt, 
+                          mytgs0[ista].ig, 
+                          fros + isrc*total_number_of_stations, 
+                          lros + isrc*total_number_of_stations, 
+                          temp_sts_data, 
+                          &istart, 
+                          &istop, 
+                          muxData);
+
+            /* weight appropriately and add */
+            for(k = 0; k < mytgs0[ista].nt; k++)
+            {
+                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
+                {
+                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
+                }
+                else
+                {
+                    sts_data[i][k] = NODATA;
+                }
+            }
+        }
+    }
+
+    free(muxData);
+    free(temp_sts_data);
+    free(fros);
+    free(lros);
+    free(mytgs0);
+
+    return sts_data;
+}
+
+/////////////////////////////////////////////////////////////////////////
+//Python gateways
+PyObject *read_mux2(PyObject *self, PyObject *args)
+{
     /*Read in mux 2 file
 
@@ -70 +489 @@
     int numSrc;
     int verbose;
-    int nsta0;
+    int total_number_of_stations;
     int number_of_selected_stations;    
     int nt;
@@ -180 +599 @@
                
     // Allocate space for return vector
-    nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
+    total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
     dt = *(double*)(file_params->data + 1*file_params->strides[0]);
     nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
@@ -286 +705 @@
 }
 
-
-float** _read_mux2(int numSrc, 
-                   char **muxFileNameArray, 
-                   float *weights, 
-                   double *params, 
-                   int *number_of_stations,
-                   long *permutation,
-                   int verbose)
-{
-    FILE *fp;
-    int nsta0, i, isrc, ista, k;
-    struct tgsrwg* mytgs0=NULL;
-    char *muxFileName;
-    int istart, istop;
-    int *fros=0, *lros=0;
-    int number_of_selected_stations;
-    float *muxData=NULL; // Suppress warning
-    long numData;
-
-    int len_sts_data;
-    float **sts_data;
-    float *temp_sts_data;
-
-    long int offset;
-
-    _read_mux2_headers(numSrc, 
-                       muxFileNameArray, 
-                       params,
-                       &fros, 
-                       &lros, 
-                       &mytgs0,
-                       &nsta0,
-                       verbose);
-    // Apply rule that an empty permutation file means 'take all stations'
-    // We could change this later by passing in None instead of the empty 
-    // permutation.
-    number_of_selected_stations = *number_of_stations;  
-    if (number_of_selected_stations == 0)
-    {
-        number_of_selected_stations = nsta0;  
-    
-        // Return possibly updated number of stations
-        *number_of_stations = nsta0;     
-      
-        // Create the Identity permutation vector
-        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
-        for (i = 0; i < number_of_selected_stations; i++)
-        {
-            permutation[i] = (long) i;  
-        }
-    }
-    
-    /*
-    printf("number_of_selected_stations = %d\n", number_of_selected_stations);
-    for (i = 0; i < number_of_selected_stations; i++) {
-      printf("permutation[%d] = %d\n", i, (int) permutation[i]);
-    }    
-    */
-    
-        
-    
-    // Make array(s) to hold demuxed data for stations given in the 
-    // permutation file 
-    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
-    if (sts_data == NULL)
-    {
-        printf("ERROR: Memory for sts_data could not be allocated.\n");
-        return NULL;
-    }
-
-    // For each selected station, allocate space for its data
-
-    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
-    for (i = 0; i < number_of_selected_stations; i++)
-    {
-        // Initialise sts_data to zero
-        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
-        if (sts_data[i] == NULL)
-        {
-            printf("ERROR: Memory for sts_data could not be allocated.\n");
-            return NULL;
-        }
-
-        ista = (int) permutation[i]; // Get global index into mux data
-    
-        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
-        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
-        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
-        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
-        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
-    }
-
-    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
-
-    /* Loop over all sources */
-    //FIXME: remove istart and istop they are not used.
-    istart = -1;
-    istop = -1;
-    for (isrc = 0; isrc < numSrc; isrc++)
-    {
-        /* Read in data block from mux2 file */
-        muxFileName = muxFileNameArray[isrc];
-        if((fp = fopen(muxFileName, "r")) == NULL)
-        {
-            fprintf(stderr, "cannot open file %s\n", muxFileName);
-            return NULL;                    
-        }
-
-        if (verbose){
-            printf("Reading mux file %s\n", muxFileName);
-        }
-
-        offset = sizeof(int) + nsta0*(sizeof(struct tgsrwg) + 2*sizeof(int));
-        fseek(fp, offset, 0);
-
-        numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
-        muxData = (float*)malloc(numData*sizeof(float));
-        fread(muxData, numData*sizeof(float), 1, fp); 
-        fclose(fp);
-
-        // loop over stations present in the permutation array 
-        //     use ista with mux data
-        //     use i with the processed data to be returned         
-        
-        for(i = 0; i < number_of_selected_stations; i++)
-        {               
-    
-            ista = (int) permutation[i]; // Get global index into mux data  
-        
-            /* fill the data0 array from the mux file, and weight it */
-            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(k = 0; k < mytgs0[ista].nt; k++)
-            {
-                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
-                {
-                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
-                }
-                else
-                {
-                    sts_data[i][k] = NODATA;
-                }
-            }
-        }
-        free(muxData);
-    }
-
-    free(temp_sts_data);
-    free(fros);
-    free(lros);
-    free(mytgs0);
-
-    return sts_data;
-}
-
-void _read_mux2_headers(int numSrc, 
-                        char **muxFileNameArray, 
-                        double *params,
-                        int** fros, 
-                        int** lros, 
-                        struct tgsrwg ** mytgs0,
-                        int* nsta0,
-                        int verbose)
-{
-    FILE *fp;
-    int nsta, isrc, ista;
-    struct tgsrwg *mytgs=0;
-    char *muxFileName;                                                                  
-    char susMuxFileName;
-    long numData;
-
-    /* Allocate space for the names and the weights and pointers to the data*/
-
-    /* Check that the input files have mux2 extension*/
-    susMuxFileName = 0;
-    for(isrc = 0; isrc < numSrc; isrc++)
-    { 
-        muxFileName = muxFileNameArray[isrc];
-        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4, 
-                     "mux2") != 0)
-        {
-            susMuxFileName = 1;
-            break;
-        }
-    }
-
-    if(susMuxFileName)
-    {
-        printf("\n**************************************************************************\n");
-        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n"); 
-        printf("   At least one input file name does not end with mux2\n");
-        printf("   Check your results carefully!\n");
-        printf("**************************************************************************\n\n");
-    }   
-
-    if (verbose)
-    {
-        printf("Reading mux header information\n");
-    }
-
-    /* Loop over all sources, read headers and check compatibility */
-    for (isrc = 0; isrc < numSrc; isrc++)
-    {
-        muxFileName = muxFileNameArray[isrc];
-
-        /* open the mux file */
-        if((fp = fopen(muxFileName, "r")) == NULL)
-        {
-            fprintf(stderr, "cannot open file %s\n", muxFileName);
-            exit(-1);  
-        }
-        
-        if (!isrc)
-        {
-            fread(nsta0, sizeof(int), 1, fp);
-        
-            *fros = (int*)malloc(*nsta0*numSrc*sizeof(int)); 
-            *lros = (int*)malloc(*nsta0*numSrc*sizeof(int));
-      
-            *mytgs0 = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
-            mytgs = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
-
-            fread(*mytgs0, *nsta0*sizeof(struct tgsrwg), 1, fp);
-        }
-        else
-        {
-            /* check that the mux files are compatible */      
-            fread(&nsta, sizeof(int), 1, fp);
-            if(nsta != *nsta0)
-            {
-                fprintf(stderr,"%s has different number of stations to %s\n", 
-                muxFileName, 
-                muxFileNameArray[0]);
-                fclose(fp);
-                exit(-1);   
-            }
-
-            fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp); 
-            
-            for (ista = 0; ista < nsta; ista++)
-            {
-                if (mytgs[ista].dt != (*mytgs0)[ista].dt)
-                {
-                    fprintf(stderr,"%s has different sampling rate to %s\n", 
-                    muxFileName, 
-                    muxFileNameArray[0]);
-                    fclose(fp);
-                    exit(-1);            
-                }   
-                if (mytgs[ista].nt != (*mytgs0)[ista].nt)
-                {
-                    fprintf(stderr,"%s has different series length to %s\n", 
-                    muxFileName, 
-                    muxFileNameArray[0]);
-                    fclose(fp);
-                    exit(-1);            
-                }
-            }
-        }
-
-        /* Read the start and stop times for this source */
-        fread(*fros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
-        fread(*lros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
-
-        /* Compute the size of the data block for this source */
-        numData = getNumData(*fros + isrc*(*nsta0), *lros + isrc*(*nsta0), (*nsta0));
-
-        /* Sanity check */
-        if (numData < 0)
-        {
-            fprintf(stderr,"Size of data block appears to be negative!\n");
-            exit(-1);        
-        }
-
-        fclose(fp);          
-    }
-
-    params[0] = (double)*nsta0;
-    params[1] = (double)(*mytgs0)[0].dt;
-    params[2] = (double)(*mytgs0)[0].nt;
-
-        free(mytgs);
-}
-
-
-/* thomas */
-void fillDataArray(int ista, int nsta, int nt, int ig, int *nst, 
-                   int *nft, float *data, int *istart_p, 
-           int *istop_p, float *muxData)
-{
-    int it, last_it, jsta;
-    long int offset=0;
-
-
-    last_it = -1;
-    /* make arrays of starting and finishing time steps for the tide gauges */
-    /* and fill them from the file */
-
-    /* update start and stop timesteps for this gauge */
-    if (nst[ista]!= -1)
-    {
-        if(*istart_p == -1)
-        {
-            *istart_p = nst[ista];
-        }
-        else
-        {
-            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
-        }
-    }
-    if (nft[ista] != -1)
-    {
-        if (*istop_p == -1)
-        {
-            *istop_p = nft[ista];
-        }
-        else
-        {
-            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
-        }
-    }     
-    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
-    {
-        /* gauge never started recording, or was outside of all grids, 
-    fill array with 0 */
-        for(it = 0; it < nt; it++)
-        {
-            data[it] = 0.0;
-        }
-    }   
-    else
-    {
-        for(it = 0; it < nt; it++)
-        {
-            last_it = it;
-            /* skip t record of data block */
-            offset++;
-            /* skip records from earlier tide gauges */
-            for(jsta = 0; jsta < ista; jsta++)
-                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
-                    offset++;
-
-            /* deal with the tide gauge at hand */
-            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
-                /* gauge is recording at this time */
-            {
-                memcpy(data + it, muxData + offset, sizeof(float));
-                offset++;
-            }
-            else if (it + 1 < nst[ista])
-            {
-                /* gauge has not yet started recording */
-                data[it] = 0.0;
-            }   
-            else
-                /* gauge has finished recording */
-            {
-                data[it] = NODATA;
-                break;
-            }
-
-            /* skip records from later tide gauges */
-            for(jsta = ista + 1; jsta < nsta; jsta++)
-                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
-                    offset++;
-        }
-
-        if(last_it < nt - 1)
-            /* the loop was exited early because the gauge had 
-        finished recording */
-            for(it = last_it+1; it < nt; it++)
-                data[it] = NODATA;
-    }
-} 
-
-char isdata(float x)
-{
-    //char value;
-    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
-    {
-       return 0;
-    }
-    else
-    {
-        return 1;  
-    }
-}
-
-
-long getNumData(const int *fros, const int *lros, const int nsta)
-/* calculates the number of data in the data block of a mux file */
-/* based on the first and last recorded output steps for each gauge */ 
-{
-    int ista, last_output_step;
-    long numData = 0;
-
-    last_output_step = 0;   
-    for(ista = 0; ista < nsta; ista++)
-        if(*(fros + ista) != -1)
-        {
-            numData += *(lros + ista) - *(fros + ista) + 1;
-            last_output_step = (last_output_step < *(lros+ista) ? 
-                            *(lros+ista):last_output_step);
-        }   
-        numData += last_output_step*nsta; /* these are the t records */
-        return numData;
-}   
-
-
-
 //-------------------------------
 // Method table for python module