Changeset 5611 for anuga_core/source/anuga/shallow_water/urs_ext.c

Timestamp:

Aug 5, 2008, 2:26:34 PM (15 years ago)

Author:

ole

Message:

Fixed urs2sts such that start and stop times are updated across mux sources

File:

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

anuga_core/source/anuga/shallow_water/urs_ext.c

@@ -57 +57 @@
         }
     }
+    
     if (nft[ista] != -1)
     {
@@ -68 +69 @@
         }
     }     
+    
     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 */
+        fill array with 0 */
         for(it = 0; it < nt; it++)
         {
@@ -91 +93 @@
             /* 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++;
@@ -116 +118 @@
         if(last_it < nt - 1)
             /* the loop was exited early because the gauge had 
-        finished recording */
+            finished recording */
             for(it = last_it+1; it < nt; it++)
                 data[it] = NODATA;
@@ -267 +269 @@
 
         /* 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);
+        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 */
@@ -291 +295 @@
 
     // Store time resolution and number of timesteps    
-    // These are the same for all stations, so 
+    // These are the same for all stations as tested above, so 
     // we take the first one.
     *delta_t = (double)mytgs0[0].dt;
@@ -312 +316 @@
     FILE *fp;
     int total_number_of_stations, i, isrc, ista, k;
-    //struct tgsrwg* mytgs0=NULL;
     char *muxFileName;
-    int istart, istop;
+    int istart=-1, istop=-1;
     int number_of_selected_stations;
     float *muxData=NULL; // Suppress warning
@@ -325 +328 @@
     long int offset;
 
-    int number_of_time_steps;
+    int number_of_time_steps, N;
     double delta_t;
-    //long numDataMax;
+    
+    // Shorthands pointing to memory blocks for each source
+    int *fros_per_source=NULL;     
+    int *lros_per_source=NULL;         
     
     _read_mux2_headers(numSrc, 
@@ -370 +376 @@
 
     // For each selected station, allocate space for its data
-
-    len_sts_data = mytgs0[0].nt + POFFSET; // Max length of each timeseries?
+    len_sts_data = number_of_time_steps + POFFSET; // Max length of each timeseries?
     for (i = 0; i < number_of_selected_stations; i++)
     {
@@ -381 +386 @@
             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];
     }
 
@@ -394 +391 @@
 
     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++)
     {
+    
+        // Shorthands to local memory
+        fros_per_source = (int*) fros + isrc*total_number_of_stations; 
+        lros_per_source = (int*) lros + isrc*total_number_of_stations;      
+            
+    
+    
         /* Read in data block from mux2 file */
         muxFileName = muxFileNameArray[isrc];
@@ -415 +417 @@
         fseek(fp, offset, 0);
 
-        numData = getNumData(fros + isrc*total_number_of_stations, 
-                             lros + isrc*total_number_of_stations, 
+        numData = getNumData(fros_per_source, 
+                             lros_per_source, 
                              total_number_of_stations);
                              
@@ -431 +433 @@
             ista = (int) permutation[i]; // Get global index into mux data  
         
-            /* fill the data0 array from the mux file, and weight it */
+            // fill the data0 array from the mux file, and weight it
             fillDataArray(ista, 
                           total_number_of_stations, 
-                          mytgs0[ista].nt, 
+                          number_of_time_steps,
                           mytgs0[ista].ig, 
-                          fros + isrc*total_number_of_stations, 
-                          lros + isrc*total_number_of_stations, 
+                          fros_per_source, 
+                          lros_per_source, 
                           temp_sts_data, 
                           &istart, 
@@ -443 +445 @@
                           muxData);
 
-            /* weight appropriately and add */
+            // Weight appropriately and add
             for(k = 0; k < mytgs0[ista].nt; k++)
             {
@@ -455 +457 @@
                 }
             }
+            
+            // Update metadata (e.g. start time and end time)
+            N = number_of_time_steps;
+            
+            
+            /*
+            printf("station =  %d, source = %d, fros = %f, lros = %f\n", 
+                   ista, 
+                   isrc, 
+                   (float) fros_per_source[ista],
+                   (float) lros_per_source[ista]);                 
+            */
+                   
+                   
+                    
+            if (isrc == 0) {
+                // Assign values for first source
+                sts_data[i][N] = (float)mytgs0[ista].geolat;
+                sts_data[i][N+1] = (float)mytgs0[ista].geolon;
+                sts_data[i][N+2] = (float)mytgs0[ista].z;
+                sts_data[i][N+3] = (float)fros_per_source[ista];
+                sts_data[i][N+4] = (float)lros_per_source[ista];
+            } else {
+                // Update first and last timesteps for subsequent sources
+                if (sts_data[i][N+3] > (float)fros_per_source[ista]) {          
+                    if (verbose) {
+                        printf("Adjusting start time for station %d and source %d",
+                               ista, isrc);
+                        printf(" from %f to %f\n", sts_data[i][N+3], (float)fros_per_source[ista]);  
+                    }
+                    sts_data[i][N+3] = (float)fros_per_source[ista];
+                }
+                
+                
+                if (sts_data[i][N+4] < (float)lros_per_source[ista]) {          
+                    if (verbose) {
+                        printf("Adjusting end time for station %d and source %d",
+                               ista, isrc);
+                        printf(" from %f to %f\n", sts_data[i][N+4], (float)lros_per_source[ista]);  
+                    }
+                    sts_data[i][N+4] = (float)lros_per_source[ista];
+                }               
+            }
         }
     }