Changeset 5560


Ignore:
Timestamp:
Jul 23, 2008, 1:38:52 PM (11 years ago)
Author:
nariman
Message:

Creation of static work arrays. Cleanup of code.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/urs_ext.c

    r5552 r5560  
    2525#define POFFSET 5 //Number of site_params
    2626
    27 void fillDataArray(int, int, int, int, int *, int *, float *,
    28            int *, int *, float *);
    29 long getNumData(const int *fros, const int *lros, const int nsta);
    30 char isdata(float);
    31 void _read_mux2_headers(int numSrc,
    32                            char **muxFileNameArray,
    33                    double *params,
    34                                    int** fros, int** lros, struct tgsrwg ** mytgs0,
    35                    int* nsta0,
    36                    int verbose);
     27static int *fros=NULL;
     28static int *lros=NULL;
     29static struct tgsrwg* mytgs0=NULL;
     30
     31static long numDataMax=0;
     32
     33/////////////////////////////////////////////////////////////////////////
     34//Auxiliary functions
     35void fillDataArray(int ista, int total_number_of_stations, int nt, int ig, int *nst,
     36                   int *nft, float *data, int *istart_p,
     37           int *istop_p, float *muxData)
     38{
     39    int it, last_it, jsta;
     40    long int offset=0;
     41
     42
     43    last_it = -1;
     44    /* make arrays of starting and finishing time steps for the tide gauges */
     45    /* and fill them from the file */
     46
     47    /* update start and stop timesteps for this gauge */
     48    if (nst[ista]!= -1)
     49    {
     50        if(*istart_p == -1)
     51        {
     52            *istart_p = nst[ista];
     53        }
     54        else
     55        {
     56            *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
     57        }
     58    }
     59    if (nft[ista] != -1)
     60    {
     61        if (*istop_p == -1)
     62        {
     63            *istop_p = nft[ista];
     64        }
     65        else
     66        {
     67            *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
     68        }
     69    }     
     70    if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
     71    {
     72        /* gauge never started recording, or was outside of all grids,
     73    fill array with 0 */
     74        for(it = 0; it < nt; it++)
     75        {
     76            data[it] = 0.0;
     77        }
     78    }   
     79    else
     80    {
     81        for(it = 0; it < nt; it++)
     82        {
     83            last_it = it;
     84            /* skip t record of data block */
     85            offset++;
     86            /* skip records from earlier tide gauges */
     87            for(jsta = 0; jsta < ista; jsta++)
     88                if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
     89                    offset++;
     90
     91            /* deal with the tide gauge at hand */
     92            if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
     93                /* gauge is recording at this time */
     94            {
     95                memcpy(data + it, muxData + offset, sizeof(float));
     96                offset++;
     97            }
     98            else if (it + 1 < nst[ista])
     99            {
     100                /* gauge has not yet started recording */
     101                data[it] = 0.0;
     102            }   
     103            else
     104                /* gauge has finished recording */
     105            {
     106                data[it] = NODATA;
     107                break;
     108            }
     109
     110            /* skip records from later tide gauges */
     111            for(jsta = ista + 1; jsta < total_number_of_stations; jsta++)
     112                if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
     113                    offset++;
     114        }
     115
     116        if(last_it < nt - 1)
     117            /* the loop was exited early because the gauge had
     118        finished recording */
     119            for(it = last_it+1; it < nt; it++)
     120                data[it] = NODATA;
     121    }
     122}
     123
     124
     125char isdata(float x)
     126{
     127    //char value;
     128    if(x < NODATA + EPSILON && NODATA < x + EPSILON)
     129    {
     130       return 0;
     131    }
     132    else
     133    {
     134        return 1; 
     135    }
     136}
     137
     138
     139long getNumData(const int *fros, const int *lros, const int total_number_of_stations)
     140/* calculates the number of data in the data block of a mux file */
     141/* based on the first and last recorded output steps for each gauge */
     142{
     143    int ista, last_output_step;
     144    long numData = 0;
     145
     146    last_output_step = 0;   
     147    for(ista = 0; ista < total_number_of_stations; ista++)
     148        if(*(fros + ista) != -1)
     149        {
     150            numData += *(lros + ista) - *(fros + ista) + 1;
     151            last_output_step = (last_output_step < *(lros+ista) ?
     152                            *(lros+ista):last_output_step);
     153        }   
     154        numData += last_output_step*total_number_of_stations; /* these are the t records */
     155        return numData;
     156}
     157
     158/////////////////////////////////////////////////////////////////////////
     159//Internal Functions
     160int _read_mux2_headers(int numSrc,
     161                        char **muxFileNameArray,
     162                        int* total_number_of_stations,
     163                                                int* number_of_time_steps,
     164                                                double* delta_t,
     165                        //long* numDataMax,
     166                        int verbose)
     167{
     168    FILE *fp;
     169    int numsta, i, j;
     170    struct tgsrwg *mytgs=0;
     171    char *muxFileName;                                                                 
     172    char susMuxFileName;
     173    long numData;
     174
     175    /* Allocate space for the names and the weights and pointers to the data*/
     176
     177    /* Check that the input files have mux2 extension*/
     178    susMuxFileName = 0;
     179    for(i = 0; i < numSrc; i++)
     180    {
     181        muxFileName = muxFileNameArray[i];
     182        if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4,
     183                     "mux2") != 0)
     184        {
     185            susMuxFileName = 1;
     186            break;
     187        }
     188    }
     189
     190    if(susMuxFileName)
     191    {
     192        printf("\n**************************************************************************\n");
     193        printf("   WARNING: This program operates only on multiplexed files in mux2 format\n");
     194        printf("   At least one input file name does not end with mux2\n");
     195        printf("   Check your results carefully!\n");
     196        printf("**************************************************************************\n\n");
     197    }   
     198
     199    if (verbose)
     200    {
     201        printf("Reading mux header information\n");
     202    }
     203
     204    /* Loop over all sources, read headers and check compatibility */
     205    for (i = 0; i < numSrc; i++)
     206    {
     207        muxFileName = muxFileNameArray[i];
     208
     209        /* open the mux file */
     210        if((fp = fopen(muxFileName, "r")) == NULL)
     211        {
     212            fprintf(stderr, "cannot open file %s\n", muxFileName);
     213            return 0; 
     214        }
     215       
     216        if (!i)
     217        {
     218            fread(total_number_of_stations, sizeof(int), 1, fp);
     219       
     220            fros = (int*)malloc(*total_number_of_stations*numSrc*sizeof(int));
     221            lros = (int*)malloc(*total_number_of_stations*numSrc*sizeof(int));
     222     
     223            mytgs0 = (struct tgsrwg*)malloc(*total_number_of_stations*sizeof(struct tgsrwg));
     224            mytgs = (struct tgsrwg*)malloc(*total_number_of_stations*sizeof(struct tgsrwg));
     225
     226            fread(mytgs0, *total_number_of_stations*sizeof(struct tgsrwg), 1, fp);
     227        }
     228        else
     229        {
     230            /* check that the mux files are compatible */     
     231            fread(&numsta, sizeof(int), 1, fp);
     232            if(numsta != *total_number_of_stations)
     233            {
     234                fprintf(stderr,"%s has different number of stations to %s\n",
     235                muxFileName,
     236                muxFileNameArray[0]);
     237                fclose(fp);
     238                return 0;   
     239            }
     240
     241            fread(mytgs, numsta*sizeof(struct tgsrwg), 1, fp);
     242           
     243            for (j = 0; j < numsta; j++)
     244            {
     245                if (mytgs[j].dt != mytgs0[j].dt)
     246                {
     247                    fprintf(stderr,"%s has different sampling rate to %s\n",
     248                    muxFileName,
     249                    muxFileNameArray[0]);
     250                    fclose(fp);
     251                    return 0;           
     252                }   
     253                if (mytgs[j].nt != mytgs0[j].nt)
     254                {
     255                    fprintf(stderr,"%s has different series length to %s\n",
     256                    muxFileName,
     257                    muxFileNameArray[0]);
     258                    fclose(fp);
     259                    return 0;           
     260                }
     261
     262                if (mytgs[j].nt != mytgs0[0].nt)
     263                {
     264                    printf("Station 0 has different series length to Station %d\n", j);
     265                }
     266            }
     267        }
     268
     269        /* Read the start and stop times for this source */
     270        fread(fros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
     271        fread(lros + i*(*total_number_of_stations), *total_number_of_stations*sizeof(int), 1, fp);
     272
     273        /* Compute the size of the data block for this source */
     274        numData = getNumData(fros + i*(*total_number_of_stations), lros + i*(*total_number_of_stations), (*total_number_of_stations));
     275
     276        /* Sanity check */
     277        if (numData < 0)
     278        {
     279            fprintf(stderr,"Size of data block appears to be negative!\n");
     280            return 0;       
     281        }
     282
     283        if (numDataMax < numData)
     284        {
     285            numDataMax = numData;
     286        }
     287
     288        fclose(fp);         
     289    }
     290
     291    *delta_t = (double)mytgs0[0].dt;
     292    *number_of_time_steps = mytgs0[0].nt;
     293
     294        free(mytgs);
     295
     296    return 1;
     297}
     298
    37299
    38300float** _read_mux2(int numSrc,
    39            char **muxFileNameArray,
    40            float *weights,
    41            double *params,
    42            int *number_of_stations,
    43            long *permutation,
    44            int verbose);
    45 
    46     PyObject *read_mux2(PyObject *self, PyObject *args){
     301                   char **muxFileNameArray,
     302                   float *weights,
     303                   double *params,
     304                   int *number_of_stations,
     305                   long *permutation,
     306                   int verbose)
     307{
     308    FILE *fp;
     309    int total_number_of_stations, i, isrc, ista, k;
     310    //struct tgsrwg* mytgs0=NULL;
     311    char *muxFileName;
     312    int istart, istop;
     313    int number_of_selected_stations;
     314    float *muxData=NULL; // Suppress warning
     315    long numData;
     316
     317    int len_sts_data;
     318    float **sts_data;
     319    float *temp_sts_data;
     320
     321    long int offset;
     322
     323    int number_of_time_steps;
     324    double delta_t;
     325    //long numDataMax;
     326   
     327    _read_mux2_headers(numSrc,
     328                       muxFileNameArray,
     329                       &total_number_of_stations,
     330                       &number_of_time_steps,
     331                       &delta_t,
     332                       verbose);
     333
     334    //FIXME: The params can be removed after the python interface is modified.
     335    params[0] = (double)total_number_of_stations;
     336    params[1] = (double)delta_t;
     337    params[2] = (double)number_of_time_steps;
     338   
     339    // Apply rule that an empty permutation file means 'take all stations'
     340    // We could change this later by passing in None instead of the empty
     341    // permutation.
     342    number_of_selected_stations = *number_of_stations; 
     343    if (number_of_selected_stations == 0)
     344    {
     345        number_of_selected_stations = total_number_of_stations; 
     346   
     347        // Return possibly updated number of stations
     348        *number_of_stations = total_number_of_stations;     
     349     
     350        // Create the Identity permutation vector
     351        permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
     352        for (i = 0; i < number_of_selected_stations; i++)
     353        {
     354            permutation[i] = (long) i; 
     355        }
     356    }
     357   
     358    // Make array(s) to hold demuxed data for stations given in the
     359    // permutation file
     360    sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
     361    if (sts_data == NULL)
     362    {
     363        printf("ERROR: Memory for sts_data could not be allocated.\n");
     364        return NULL;
     365    }
     366
     367    // For each selected station, allocate space for its data
     368
     369    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
     370    for (i = 0; i < number_of_selected_stations; i++)
     371    {
     372        // Initialise sts_data to zero
     373        sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
     374        if (sts_data[i] == NULL)
     375        {
     376            printf("ERROR: Memory for sts_data could not be allocated.\n");
     377            return NULL;
     378        }
     379
     380        ista = (int) permutation[i]; // Get global index into mux data
     381   
     382        sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
     383        sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
     384        sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
     385        sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
     386        sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
     387    }
     388
     389    temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
     390
     391    muxData = (float*)malloc(numDataMax*sizeof(float));
     392    /* Loop over all sources */
     393    //FIXME: remove istart and istop they are not used.
     394    istart = -1;
     395    istop = -1;
     396    for (isrc = 0; isrc < numSrc; isrc++)
     397    {
     398        /* Read in data block from mux2 file */
     399        muxFileName = muxFileNameArray[isrc];
     400        if((fp = fopen(muxFileName, "r")) == NULL)
     401        {
     402            fprintf(stderr, "cannot open file %s\n", muxFileName);
     403            return NULL;                   
     404        }
     405
     406        if (verbose){
     407            printf("Reading mux file %s\n", muxFileName);
     408        }
     409
     410        offset = sizeof(int) + total_number_of_stations*(sizeof(struct tgsrwg) + 2*sizeof(int));
     411        fseek(fp, offset, 0);
     412
     413        numData = getNumData(fros + isrc*total_number_of_stations, lros + isrc*total_number_of_stations, total_number_of_stations);
     414        fread(muxData, numData*sizeof(float), 1, fp);
     415        fclose(fp);
     416
     417        // loop over stations present in the permutation array
     418        //     use ista with mux data
     419        //     use i with the processed data to be returned         
     420       
     421        for(i = 0; i < number_of_selected_stations; i++)
     422        {               
     423   
     424            ista = (int) permutation[i]; // Get global index into mux data 
     425       
     426            /* fill the data0 array from the mux file, and weight it */
     427            fillDataArray(ista,
     428                          total_number_of_stations,
     429                          mytgs0[ista].nt,
     430                          mytgs0[ista].ig,
     431                          fros + isrc*total_number_of_stations,
     432                          lros + isrc*total_number_of_stations,
     433                          temp_sts_data,
     434                          &istart,
     435                          &istop,
     436                          muxData);
     437
     438            /* weight appropriately and add */
     439            for(k = 0; k < mytgs0[ista].nt; k++)
     440            {
     441                if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
     442                {
     443                    sts_data[i][k] += temp_sts_data[k] * weights[isrc];
     444                }
     445                else
     446                {
     447                    sts_data[i][k] = NODATA;
     448                }
     449            }
     450        }
     451    }
     452
     453    free(muxData);
     454    free(temp_sts_data);
     455    free(fros);
     456    free(lros);
     457    free(mytgs0);
     458
     459    return sts_data;
     460}
     461
     462/////////////////////////////////////////////////////////////////////////
     463//Python gateways
     464PyObject *read_mux2(PyObject *self, PyObject *args)
     465{
    47466    /*Read in mux 2 file
    48467
     
    70489    int numSrc;
    71490    int verbose;
    72     int nsta0;
     491    int total_number_of_stations;
    73492    int number_of_selected_stations;   
    74493    int nt;
     
    180599               
    181600    // Allocate space for return vector
    182     nsta0 = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
     601    total_number_of_stations = (int)*(double*)(file_params->data + 0*file_params->strides[0]);
    183602    dt = *(double*)(file_params->data + 1*file_params->strides[0]);
    184603    nt = (int)*(double*)(file_params->data + 2*file_params->strides[0]);
     
    286705}
    287706
    288 
    289 float** _read_mux2(int numSrc,
    290                    char **muxFileNameArray,
    291                    float *weights,
    292                    double *params,
    293                    int *number_of_stations,
    294                    long *permutation,
    295                    int verbose)
    296 {
    297     FILE *fp;
    298     int nsta0, i, isrc, ista, k;
    299     struct tgsrwg* mytgs0=NULL;
    300     char *muxFileName;
    301     int istart, istop;
    302     int *fros=0, *lros=0;
    303     int number_of_selected_stations;
    304     float *muxData=NULL; // Suppress warning
    305     long numData;
    306 
    307     int len_sts_data;
    308     float **sts_data;
    309     float *temp_sts_data;
    310 
    311     long int offset;
    312 
    313     _read_mux2_headers(numSrc,
    314                        muxFileNameArray,
    315                        params,
    316                        &fros,
    317                        &lros,
    318                        &mytgs0,
    319                        &nsta0,
    320                        verbose);
    321     // Apply rule that an empty permutation file means 'take all stations'
    322     // We could change this later by passing in None instead of the empty
    323     // permutation.
    324     number_of_selected_stations = *number_of_stations; 
    325     if (number_of_selected_stations == 0)
    326     {
    327         number_of_selected_stations = nsta0; 
    328    
    329         // Return possibly updated number of stations
    330         *number_of_stations = nsta0;     
    331      
    332         // Create the Identity permutation vector
    333         permutation = (long *) malloc(number_of_selected_stations*sizeof(long));
    334         for (i = 0; i < number_of_selected_stations; i++)
    335         {
    336             permutation[i] = (long) i; 
    337         }
    338     }
    339    
    340     /*
    341     printf("number_of_selected_stations = %d\n", number_of_selected_stations);
    342     for (i = 0; i < number_of_selected_stations; i++) {
    343       printf("permutation[%d] = %d\n", i, (int) permutation[i]);
    344     }   
    345     */
    346    
    347        
    348    
    349     // Make array(s) to hold demuxed data for stations given in the
    350     // permutation file
    351     sts_data = (float**)malloc(number_of_selected_stations*sizeof(float*));
    352     if (sts_data == NULL)
    353     {
    354         printf("ERROR: Memory for sts_data could not be allocated.\n");
    355         return NULL;
    356     }
    357 
    358     // For each selected station, allocate space for its data
    359 
    360     len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
    361     for (i = 0; i < number_of_selected_stations; i++)
    362     {
    363         // Initialise sts_data to zero
    364         sts_data[i] = (float*)calloc(len_sts_data, sizeof(float));
    365         if (sts_data[i] == NULL)
    366         {
    367             printf("ERROR: Memory for sts_data could not be allocated.\n");
    368             return NULL;
    369         }
    370 
    371         ista = (int) permutation[i]; // Get global index into mux data
    372    
    373         sts_data[i][mytgs0[0].nt] = (float)mytgs0[ista].geolat;
    374         sts_data[i][mytgs0[0].nt + 1] = (float)mytgs0[ista].geolon;
    375         sts_data[i][mytgs0[0].nt + 2] = (float)mytgs0[ista].z;
    376         sts_data[i][mytgs0[0].nt + 3] = (float)fros[ista];
    377         sts_data[i][mytgs0[0].nt + 4] = (float)lros[ista];
    378     }
    379 
    380     temp_sts_data = (float*)calloc(len_sts_data, sizeof(float));
    381 
    382     /* Loop over all sources */
    383     //FIXME: remove istart and istop they are not used.
    384     istart = -1;
    385     istop = -1;
    386     for (isrc = 0; isrc < numSrc; isrc++)
    387     {
    388         /* Read in data block from mux2 file */
    389         muxFileName = muxFileNameArray[isrc];
    390         if((fp = fopen(muxFileName, "r")) == NULL)
    391         {
    392             fprintf(stderr, "cannot open file %s\n", muxFileName);
    393             return NULL;                   
    394         }
    395 
    396         if (verbose){
    397             printf("Reading mux file %s\n", muxFileName);
    398         }
    399 
    400         offset = sizeof(int) + nsta0*(sizeof(struct tgsrwg) + 2*sizeof(int));
    401         fseek(fp, offset, 0);
    402 
    403         numData = getNumData(fros + isrc*nsta0, lros + isrc*nsta0, nsta0);
    404         muxData = (float*)malloc(numData*sizeof(float));
    405         fread(muxData, numData*sizeof(float), 1, fp);
    406         fclose(fp);
    407 
    408         // loop over stations present in the permutation array
    409         //     use ista with mux data
    410         //     use i with the processed data to be returned         
    411        
    412         for(i = 0; i < number_of_selected_stations; i++)
    413         {               
    414    
    415             ista = (int) permutation[i]; // Get global index into mux data 
    416        
    417             /* fill the data0 array from the mux file, and weight it */
    418             fillDataArray(ista,
    419                           nsta0,
    420                           mytgs0[ista].nt,
    421                           mytgs0[ista].ig,
    422                           fros + isrc*nsta0,
    423                           lros + isrc*nsta0,
    424                           temp_sts_data,
    425                           &istart,
    426                           &istop,
    427                           muxData);
    428 
    429             /* weight appropriately and add */
    430             for(k = 0; k < mytgs0[ista].nt; k++)
    431             {
    432                 if((isdata(sts_data[i][k])) && isdata(temp_sts_data[k]))
    433                 {
    434                     sts_data[i][k] += temp_sts_data[k] * weights[isrc];
    435                 }
    436                 else
    437                 {
    438                     sts_data[i][k] = NODATA;
    439                 }
    440             }
    441         }
    442         free(muxData);
    443     }
    444 
    445     free(temp_sts_data);
    446     free(fros);
    447     free(lros);
    448     free(mytgs0);
    449 
    450     return sts_data;
    451 }
    452 
    453 void _read_mux2_headers(int numSrc,
    454                         char **muxFileNameArray,
    455                         double *params,
    456                         int** fros,
    457                         int** lros,
    458                         struct tgsrwg ** mytgs0,
    459                         int* nsta0,
    460                         int verbose)
    461 {
    462     FILE *fp;
    463     int nsta, isrc, ista;
    464     struct tgsrwg *mytgs=0;
    465     char *muxFileName;                                                                 
    466     char susMuxFileName;
    467     long numData;
    468 
    469     /* Allocate space for the names and the weights and pointers to the data*/
    470 
    471     /* Check that the input files have mux2 extension*/
    472     susMuxFileName = 0;
    473     for(isrc = 0; isrc < numSrc; isrc++)
    474     {
    475         muxFileName = muxFileNameArray[isrc];
    476         if(!susMuxFileName && strcmp(muxFileName + strlen(muxFileName) - 4,
    477                      "mux2") != 0)
    478         {
    479             susMuxFileName = 1;
    480             break;
    481         }
    482     }
    483 
    484     if(susMuxFileName)
    485     {
    486         printf("\n**************************************************************************\n");
    487         printf("   WARNING: This program operates only on multiplexed files in mux2 format\n");
    488         printf("   At least one input file name does not end with mux2\n");
    489         printf("   Check your results carefully!\n");
    490         printf("**************************************************************************\n\n");
    491     }   
    492 
    493     if (verbose)
    494     {
    495         printf("Reading mux header information\n");
    496     }
    497 
    498     /* Loop over all sources, read headers and check compatibility */
    499     for (isrc = 0; isrc < numSrc; isrc++)
    500     {
    501         muxFileName = muxFileNameArray[isrc];
    502 
    503         /* open the mux file */
    504         if((fp = fopen(muxFileName, "r")) == NULL)
    505         {
    506             fprintf(stderr, "cannot open file %s\n", muxFileName);
    507             exit(-1); 
    508         }
    509        
    510         if (!isrc)
    511         {
    512             fread(nsta0, sizeof(int), 1, fp);
    513        
    514             *fros = (int*)malloc(*nsta0*numSrc*sizeof(int));
    515             *lros = (int*)malloc(*nsta0*numSrc*sizeof(int));
    516      
    517             *mytgs0 = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
    518             mytgs = (struct tgsrwg*)malloc(*nsta0*sizeof(struct tgsrwg));
    519 
    520             fread(*mytgs0, *nsta0*sizeof(struct tgsrwg), 1, fp);
    521         }
    522         else
    523         {
    524             /* check that the mux files are compatible */     
    525             fread(&nsta, sizeof(int), 1, fp);
    526             if(nsta != *nsta0)
    527             {
    528                 fprintf(stderr,"%s has different number of stations to %s\n",
    529                 muxFileName,
    530                 muxFileNameArray[0]);
    531                 fclose(fp);
    532                 exit(-1);   
    533             }
    534 
    535             fread(mytgs, nsta*sizeof(struct tgsrwg), 1, fp);
    536            
    537             for (ista = 0; ista < nsta; ista++)
    538             {
    539                 if (mytgs[ista].dt != (*mytgs0)[ista].dt)
    540                 {
    541                     fprintf(stderr,"%s has different sampling rate to %s\n",
    542                     muxFileName,
    543                     muxFileNameArray[0]);
    544                     fclose(fp);
    545                     exit(-1);           
    546                 }   
    547                 if (mytgs[ista].nt != (*mytgs0)[ista].nt)
    548                 {
    549                     fprintf(stderr,"%s has different series length to %s\n",
    550                     muxFileName,
    551                     muxFileNameArray[0]);
    552                     fclose(fp);
    553                     exit(-1);           
    554                 }
    555             }
    556         }
    557 
    558         /* Read the start and stop times for this source */
    559         fread(*fros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
    560         fread(*lros + isrc*(*nsta0), *nsta0*sizeof(int), 1, fp);
    561 
    562         /* Compute the size of the data block for this source */
    563         numData = getNumData(*fros + isrc*(*nsta0), *lros + isrc*(*nsta0), (*nsta0));
    564 
    565         /* Sanity check */
    566         if (numData < 0)
    567         {
    568             fprintf(stderr,"Size of data block appears to be negative!\n");
    569             exit(-1);       
    570         }
    571 
    572         fclose(fp);         
    573     }
    574 
    575     params[0] = (double)*nsta0;
    576     params[1] = (double)(*mytgs0)[0].dt;
    577     params[2] = (double)(*mytgs0)[0].nt;
    578 
    579         free(mytgs);
    580 }
    581 
    582 
    583 /* thomas */
    584 void fillDataArray(int ista, int nsta, int nt, int ig, int *nst,
    585                    int *nft, float *data, int *istart_p,
    586            int *istop_p, float *muxData)
    587 {
    588     int it, last_it, jsta;
    589     long int offset=0;
    590 
    591 
    592     last_it = -1;
    593     /* make arrays of starting and finishing time steps for the tide gauges */
    594     /* and fill them from the file */
    595 
    596     /* update start and stop timesteps for this gauge */
    597     if (nst[ista]!= -1)
    598     {
    599         if(*istart_p == -1)
    600         {
    601             *istart_p = nst[ista];
    602         }
    603         else
    604         {
    605             *istart_p = ((nst[ista] < *istart_p) ? nst[ista] : *istart_p);
    606         }
    607     }
    608     if (nft[ista] != -1)
    609     {
    610         if (*istop_p == -1)
    611         {
    612             *istop_p = nft[ista];
    613         }
    614         else
    615         {
    616             *istop_p = ((nft[ista] < *istop_p) ? nft[ista] : *istop_p);
    617         }
    618     }     
    619     if (ig == -1 || nst[ista] == -1) /* currently ig==-1 => nst[ista]==-1 */
    620     {
    621         /* gauge never started recording, or was outside of all grids,
    622     fill array with 0 */
    623         for(it = 0; it < nt; it++)
    624         {
    625             data[it] = 0.0;
    626         }
    627     }   
    628     else
    629     {
    630         for(it = 0; it < nt; it++)
    631         {
    632             last_it = it;
    633             /* skip t record of data block */
    634             offset++;
    635             /* skip records from earlier tide gauges */
    636             for(jsta = 0; jsta < ista; jsta++)
    637                 if(it + 1 >= nst[jsta] && it + 1 <= nft[jsta])
    638                     offset++;
    639 
    640             /* deal with the tide gauge at hand */
    641             if(it + 1 >= nst[ista] && it + 1 <= nft[ista])
    642                 /* gauge is recording at this time */
    643             {
    644                 memcpy(data + it, muxData + offset, sizeof(float));
    645                 offset++;
    646             }
    647             else if (it + 1 < nst[ista])
    648             {
    649                 /* gauge has not yet started recording */
    650                 data[it] = 0.0;
    651             }   
    652             else
    653                 /* gauge has finished recording */
    654             {
    655                 data[it] = NODATA;
    656                 break;
    657             }
    658 
    659             /* skip records from later tide gauges */
    660             for(jsta = ista + 1; jsta < nsta; jsta++)
    661                 if(it + 1 >= nst[jsta] && it+1 <= nft[jsta])
    662                     offset++;
    663         }
    664 
    665         if(last_it < nt - 1)
    666             /* the loop was exited early because the gauge had
    667         finished recording */
    668             for(it = last_it+1; it < nt; it++)
    669                 data[it] = NODATA;
    670     }
    671 }
    672 
    673 char isdata(float x)
    674 {
    675     //char value;
    676     if(x < NODATA + EPSILON && NODATA < x + EPSILON)
    677     {
    678        return 0;
    679     }
    680     else
    681     {
    682         return 1; 
    683     }
    684 }
    685 
    686 
    687 long getNumData(const int *fros, const int *lros, const int nsta)
    688 /* calculates the number of data in the data block of a mux file */
    689 /* based on the first and last recorded output steps for each gauge */
    690 {
    691     int ista, last_output_step;
    692     long numData = 0;
    693 
    694     last_output_step = 0;   
    695     for(ista = 0; ista < nsta; ista++)
    696         if(*(fros + ista) != -1)
    697         {
    698             numData += *(lros + ista) - *(fros + ista) + 1;
    699             last_output_step = (last_output_step < *(lros+ista) ?
    700                             *(lros+ista):last_output_step);
    701         }   
    702         numData += last_output_step*nsta; /* these are the t records */
    703         return numData;
    704 }   
    705 
    706 
    707 
    708707//-------------------------------
    709708// Method table for python module
Note: See TracChangeset for help on using the changeset viewer.