Changeset 272


Ignore:
Timestamp:
Sep 6, 2004, 11:25:50 AM (20 years ago)
Author:
ole
Message:

Coded update as C implementation and verified improvement.
Removed oldstyle obsolete balanced limiter

Location:
inundation/ga/storm_surge/pyvolution
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/HUSK.txt

    r269 r272  
    11
     2       
     3       
    24Perhaps control old/new style with minimum_allowed height
     5In any case get rid of oldstyle.
     6Then optimise protect_against.....
    37       
    48Test friction term!

  • inundation/ga/storm_surge/pyvolution/config.py

    r268 r272  
    55
    66default_boundary_tag = 'exterior'
    7 
    8 newstyle = False #Flase meqans that we are equivalent to previous version pyvolution2
    97
    108
    … …  
    5755
    5856
    59 
    6057use_extensions = True   #Try to use C-extensions
    6158#use_extensions = False   #Do not use C-extensions
    6259
    6360use_psyco = True  #Use psyco optimisations
    64 #use_psyco = False  #Do not use psyco optimisations
     61use_psyco = False  #Do not use psyco optimisations
    6562
    6663

  • inundation/ga/storm_surge/pyvolution/domain.py

    r269 r272  
    7272
    7373        #Defaults
    74         from config import max_smallsteps, beta, newstyle, epsilon
     74        from config import max_smallsteps, beta, epsilon
    7575        self.beta = beta
    7676        self.epsilon = epsilon       
    77         self.newstyle = newstyle
    7877        self.default_order = 1
    7978        self.order = self.default_order

  • inundation/ga/storm_surge/pyvolution/quantity.py

    r265 r272  
    226226
    227227    def update(self, timestep):
    228         """Update centroid values based on values stored in
    229         explicit_update and semi_implicit_update as well as given timestep
    230         """
    231 
    232         from Numeric import sum, equal, ones, Float
    233        
    234         N = self.centroid_values.shape[0]
    235        
    236         #Explicit updates
    237         self.centroid_values += timestep*self.explicit_update
    238            
    239         #Semi implicit updates
    240         denominator = ones(N, Float)-timestep*self.semi_implicit_update
    241 
    242         if sum(equal(denominator, 0.0)) > 0.0:
    243             msg = 'Zero division in semi implicit update. Call Stephen :-)'
    244             raise msg
    245         else:
    246             #Update conserved_quantities from semi implicit updates
    247             self.centroid_values /= denominator
     228        #Call correct module function
     229        #(either from this module or C-extension)
     230        return update(self, timestep)
    248231
    249232
    … …  
    278261        extrapolate_second_order(self)
    279262
     263
     264def update(quantity, timestep):   
     265    """Update centroid values based on values stored in
     266    explicit_update and semi_implicit_update as well as given timestep
     267    """
     268   
     269    from Numeric import sum, equal, ones, Float
     270       
     271    N = quantity.centroid_values.shape[0]
     272       
     273    #Explicit updates
     274    quantity.centroid_values += timestep*quantity.explicit_update
     275           
     276    #Semi implicit updates
     277    denominator = ones(N, Float)-timestep*quantity.semi_implicit_update
     278
     279    if sum(equal(denominator, 0.0)) > 0.0:
     280        msg = 'Zero division in semi implicit update. Call Stephen :-)'
     281        raise msg
     282    else:
     283        #Update conserved_quantities from semi implicit updates
     284        quantity.centroid_values /= denominator
    280285
    281286
    … …  
    462467       
    463468    from quantity_ext import limit, compute_gradients,\
    464     extrapolate_second_order, interpolate_from_vertices_to_edges
    465 
     469    extrapolate_second_order, interpolate_from_vertices_to_edges, update
     470

  • inundation/ga/storm_surge/pyvolution/quantity_ext.c

    r265 r272  
    156156  return 0;
    157157}
     158
     159int _update(int N,
     160            double timestep,
     161            double* centroid_values,
     162            double* explicit_update,
     163            double* semi_implicit_update) {
     164  //Update centroid values based on values stored in
     165  //explicit_update and semi_implicit_update as well as given timestep
     166 
     167 
     168  int k;
     169  double denominator;
     170 
     171  //Explicit updates
     172  for (k=0; k<N; k++) {
     173    centroid_values[k] += timestep*explicit_update[k];
     174  }
     175           
     176  //Semi implicit updates
     177  for (k=0; k<N; k++) { 
     178    denominator = 1.0 - timestep*semi_implicit_update[k]; 
     179   
     180    if (denominator == 0.0) {
     181      return -1;
     182    } else {
     183      //Update conserved_quantities from semi implicit updates
     184      centroid_values[k] /= denominator;
     185    }
     186  }
     187  return 0;
     188}
    158189                     
    159190 
    160191/////////////////////////////////////////////////
    161192// Gateways to Python
     193PyObject *update(PyObject *self, PyObject *args) {
     194 
     195  PyObject *quantity;
     196  PyArrayObject *centroid_values, *explicit_update, *semi_implicit_update;
     197 
     198  double timestep; 
     199  int N, err;
     200 
     201  // Convert Python arguments to C 
     202  if (!PyArg_ParseTuple(args, "Od", &quantity, &timestep))
     203    return NULL;
     204
     205  centroid_values = (PyArrayObject*)
     206    PyObject_GetAttrString(quantity, "centroid_values");           
     207  if (!centroid_values) return NULL;               
     208
     209  explicit_update = (PyArrayObject*)
     210    PyObject_GetAttrString(quantity, "explicit_update");   
     211  if (!explicit_update) return NULL;         
     212 
     213  semi_implicit_update = (PyArrayObject*)
     214    PyObject_GetAttrString(quantity, "semi_implicit_update");   
     215  if (!semi_implicit_update) return NULL;           
     216
     217  N = centroid_values -> dimensions[0]; 
     218 
     219 
     220  err = _update(N, timestep,
     221                (double*) centroid_values -> data,
     222                (double*) explicit_update -> data,
     223                (double*) semi_implicit_update -> data);               
     224               
     225                     
     226  if (err != 0) {
     227    PyErr_SetString(PyExc_RuntimeError,
     228                    "Zero division in semi implicit update - call Stephen :)");
     229    return NULL;
     230  }                 
     231 
     232  //Release and return               
     233  Py_DECREF(centroid_values);   
     234  Py_DECREF(explicit_update);
     235  Py_DECREF(semi_implicit_update);   
     236
     237  return Py_BuildValue("");
     238}
    162239
    163240
    … …  
    471548static struct PyMethodDef MethodTable[] = {
    472549  {"limit", limit, METH_VARARGS, "Print out"},   
     550  {"update", update, METH_VARARGS, "Print out"},     
    473551  {"compute_gradients", compute_gradients, METH_VARARGS, "Print out"},     
    474552  {"extrapolate_second_order", extrapolate_second_order,

  • inundation/ga/storm_surge/pyvolution/run_profile.py

    r268 r272  
    3030domain.default_order = 2
    3131#domain.visualise = True
    32 domain.newstyle = False  #For comparison
    3332
    3433

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r271 r272  
    5757        self.z = self.quantities['elevation']           
    5858        self.eta = self.quantities['friction']                 
    59        
    60    
     59
    6160    def check_integrity(self):
    6261        Generic_domain.check_integrity(self)
    … …  
    420419    #Update
    421420    for k in range(domain.number_of_elements):
    422         if domain.newstyle:
    423             if hc[k] < domain.minimum_allowed_height:       
    424                 if hc[k] < 0.0:
    425                     #Control level and height
    426                     wc[k] = zc[k] #; hc[k] = 0.0
    427 
    428                 #Control momentum
    429                 xmomc[k] = ymomc[k] = 0.0
    430         else:       
     421        if hc[k] < domain.minimum_allowed_height:       
    431422            if hc[k] < 0.0:
    432423                #Control level and height
    433                 wc[k] = zc[k]  #; hc[k] = 0.0
    434 
    435                 #Control momentum
    436                 xmomc[k] = ymomc[k] = 0.0
    437 
     424                wc[k] = zc[k]
     425
     426            #Control momentum
     427            xmomc[k] = ymomc[k] = 0.0
    438428
    439429def extrapolate_first_order(domain):
    … …  
    461451
    462452
    463     if not domain.newstyle:
    464         #Protect against infinitesimal heights and high speeds at_centroid
    465         #FIXME: To be phased out
    466        
    467         for k in range(domain.number_of_elements):
    468             #Protect against infinitesimal heights and high velocities
    469             if hc[k] < domain.minimum_allowed_height:
    470                 #Control level and height
    471                 if hc[k] < 0.0:
    472                     wc[k] = zc[k]; hc[k] = 0.0#
    473          
    474                 #Control momentum
    475                 xmomc[k] = ymomc[k] = 0.0
    476                
    477    
    478453    #Update conserved quantities using straight first order
    479454    for name in domain.conserved_quantities:
    … …  
    483458
    484459
    485     if domain.newstyle:
    486         balance_deep_and_shallow(domain)                      #This will be better
    487     else:   
    488         old_first_order_balancing_of_deep_and_shallow(domain) #Tests will pass
     460
     461    balance_deep_and_shallow(domain)
     462
    489463
    490464   
    … …  
    562536        hmin = min( hv[k,:] )
    563537
    564         print dz, hmin
    565538       
    566539        #Create alpha in [0,1], where alpha==0 means using shallow
    … …  
    968941            return self.W(x,y) + self.h
    969942
    970 #STUFF
    971 
    972 
    973 def old_first_order_balancing_of_deep_and_shallow(domain):
    974     #FIXME: Will be obsolete, but keep comments from this one.
    975    
    976     #Computed weighted balance between constant levels and and
    977     #levels parallel to the bed elevation.
    978     wc = domain.quantities['level'].centroid_values
    979     zc = domain.quantities['elevation'].centroid_values
    980     xmomc = domain.quantities['xmomentum'].centroid_values
    981     ymomc = domain.quantities['ymomentum'].centroid_values           
    982     hc = wc - zc   #Water depths at centroids
    983    
    984     wv = domain.quantities['level'].vertex_values
    985     zv = domain.quantities['elevation'].vertex_values
    986    
    987     for k in range(domain.number_of_elements):               
    988                
    989         #Compute maximal variation in bed elevation
    990         z_range = max(abs(zv[k,0]-zc[k]),
    991                       abs(zv[k,1]-zc[k]),
    992                       abs(zv[k,2]-zc[k]))
    993 
    994 
    995         #Weighted balance between stage parallel to bed elevation
    996         #(wvi = zvi + hc) and constant stage (wvi = wc = zc+hc)
    997         #where i=0,1,2 denotes the vertex ids
    998         #
    999         #It follows that
    1000         #  wvi = (1-alpha)*(zvi+hc) + alpha*(zc+hc) =
    1001         #  (1-alpha)*zvi + alpha*zc + hc =
    1002         #  zvi + hc + alpha*(zc-zvi)
    1003         #
    1004         #where alpha in [0,1] and defined as the ratio between hc and
    1005         #the maximal difference from zc to zv0, zv1 and zv2
    1006         #
    1007         #Mathematically the following could be continued on using hc as
    1008         #  wvi =
    1009         #  zvi + hc + alpha*(zc+hc-zvi-hc) =
    1010         #  zvi + hc + alpha*(hvi-hc)
    1011         #since hvi = zc+hc-zvi in the constant case
    1012 
    1013        
    1014         if z_range > 0.0:
    1015             alpha = min(hc[k]/z_range, 1.0)
    1016         else:
    1017             alpha = 1.0
    1018            
    1019         #Update water levels
    1020         for i in range(3):
    1021             #FIXME: Use the original first-order one first, then switch
    1022             wv[k,i] = zv[k,i] + hc[k] + alpha*(zc[k]-zv[k,i])
    1023             #wv[k,i] = zv[k,i] + hc[k] + alpha*(hv[k,i]-hc[k]) #Requires updated hv!!
    1024 
    1025         #FIXME: What about alpha weighting of momentum??
    1026 
    1027943
    1028944##############################################

  • inundation/ga/storm_surge/pyvolution/show_balanced_limiters.py

    r270 r272  
    3838domain.visualise = True
    3939domain.default_order = 2
    40 domain.newstyle = True  #True is better
    4140
    4241

  • inundation/ga/storm_surge/pyvolution/test_shallow_water.py

    r240 r272  
    1212    def setUp(self):
    1313        pass
    14         #print "  Setting up"
    1514       
    1615    def tearDown(self):
    1716        pass
    18         #print "  Tearing down"
    1917
    2018    def test_rotate(self):
    … …  
    823821        domain.order = 1       
    824822        domain.distribute_to_vertices_and_edges()
    825         assert allclose(L[1], [0.19999999, 20.05, 20.05])
    826 
     823        ##assert allclose(L[1], [0.19999999, 20.05, 20.05])
     824        assert allclose(L[1], [0.1, 20.1, 20.1])
     825       
    827826        domain.order = 2       
    828827        domain.distribute_to_vertices_and_edges()
    … …  
    865864
    866865        #print E
    867         domain.order = 1       
     866        domain.order = 1
    868867        domain.distribute_to_vertices_and_edges()
    869         assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
     868        ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143])
     869        assert allclose(L[1], [4.1, 16.1, 20.1])       
    870870
    871871        domain.order = 2       
    … …  
    10851085        domain.visualise = False
    10861086        domain.default_order=2
     1087        #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance
    10871088
    10881089        # Boundary conditions
    … …  
    11061107       
    11071108        assert allclose(domain.quantities['level'].vertex_values[:4,0],
    1108                         [0.00101913,0.05352143,0.00104852,0.05354394])
     1109                        [0.00101913,0.05352143,0.00104852,0.05354394])
     1110       
    11091111        assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    1110                         [0.00090581,0.03685719,0.00088303,0.03687313])
     1112                        [ 0.00064835, 0.03685719, 0.00085073, 0.03687313])
     1113
     1114        #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     1115        #                [0.00090581,0.03685719,0.00088303,0.03687313])
     1116       
    11111117        assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
    11121118                        [-0.00139463,0.0006156,-0.00060364,0.00061827])
    … …  
    12121218            #print 'C19=', domain.quantities['xmomentum'].centroid_values[19]
    12131219
    1214             assert allclose(domain.quantities['xmomentum'].centroid_values[17],
    1215                             0.00028606084)
     1220            #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084)
     1221            ##print domain.quantities['xmomentum'].centroid_values[17], V
     1222            ##print
     1223            if not V:           
     1224                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0)
     1225            else: 
     1226                assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084)           
    12161227
    12171228            import copy
    … …  
    12211232            domain.distribute_to_vertices_and_edges()
    12221233
    1223             assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
     1234            #assert allclose(domain.quantities['xmomentum'].centroid_values, XX)
    12241235
    12251236            assert allclose(domain.quantities['xmomentum'].centroid_values[17],
    1226                             0.00028606084)
     1237                            0.0)
    12271238
    12281239
    … …  
    12321243            assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0],
    12331244                            [-0.00139463,0.0006156,-0.00060364,0.00061827])
     1245
     1246
     1247            assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     1248                            [0.00064835,0.03685719,0.00085073,0.03687313])
     1249
     1250
     1251            #NB NO longer relvant:
    12341252
    12351253            #This was the culprit. First triangles vertex 0 had an
    … …  
    12371255            #third triangle had 0.00850733 instead of 0.00088303
    12381256            #print domain.quantities['xmomentum'].vertex_values[:4,0]
    1239        
    1240             assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
    1241                             [0.00090581,0.03685719,0.00088303,0.03687313])       
     1257
     1258            #print domain.quantities['xmomentum'].vertex_values[:4,0]
     1259            #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0],
     1260            #                [0.00090581,0.03685719,0.00088303,0.03687313])       
    12421261
    12431262
    … …  
    12771296            pass# domain.write_time()
    12781297
    1279         #Data from earlier version of pyvolution
    1280         assert allclose(domain.min_timestep, 0.0423354136479)
    1281         assert allclose(domain.max_timestep, 0.0479397182296)
    1282         assert allclose(domain.quantities['level'].centroid_values,
    1283                         [0.00974056, 0.02131643, 0.01373269, 0.02139684,
    1284                          0.01373269, 0.02139684, 0.01373269, 0.02139684,
    1285                          0.01373269, 0.02139684, 0.01386079, 0.0254346,
    1286                          -0.04562579, -0.025253, -0.04154658, -0.02512009,
    1287                          -0.04154658, -0.02512009, -0.04154658, -0.02512009,
    1288                          -0.04154658, -0.02512009, -0.04140825, -0.0210463,
    1289                          -0.10118135, -0.08080855, -0.09710213, -0.08067564,
    1290                          -0.09710213, -0.08067564, -0.09710213, -0.08067564,
    1291                          -0.09710213, -0.08067564, -0.0969638 , -0.07660185,
    1292                          -0.1567369, -0.13636411, -0.15265769, -0.1362312,
    1293                          -0.15265769, -0.1362312, -0.15265769, -0.1362312,
    1294                          -0.15265769, -0.1362312, -0.15251936, -0.13215741,
    1295                          -0.21229246, -0.19191967, -0.20821325, -0.19178675,
    1296                          -0.20821325, -0.19178675, -0.20821325, -0.19178675,
    1297                          -0.20821325, -0.19178675, -0.20807491, -0.18771296,
    1298                          -0.25875433, -0.24722502, -0.25470346, -0.24709274,
    1299                          -0.25470346, -0.24709274, -0.25470346, -0.24709274,
    1300                          -0.25470346, -0.24709274, -0.25460675, -0.24309962])
    1301 
     1298        assert allclose(domain.min_timestep, 0.050010003001)
     1299        assert allclose(domain.max_timestep, 0.050010003001)
    13021300
    13031301       
    … …  
    13351333
    13361334        #Data from earlier version of pyvolution
     1335        #print domain.quantities['level'].centroid_values
     1336       
    13371337        assert allclose(domain.quantities['level'].centroid_values,
    1338                         [-0.03558044, -0.01758186, -0.03457868,
    1339                          -0.01757463, -0.03400681, -0.01765817,
    1340                          -0.0334924, -0.01763917, -0.03304351,
    1341                          -0.01745604, -0.03259134, -0.01690249,
    1342                          -0.08343238, -0.06844054, -0.08114626,
    1343                          -0.06776272, -0.07985917, -0.0670156,
    1344                          -0.07857347, -0.06636824, -0.07751765,
    1345                          -0.06571005, -0.07694892, -0.06485356,
    1346                          -0.12816762, -0.11399292, -0.12388354,
    1347                          -0.11298543, -0.12146159, -0.11160487,
    1348                          -0.11890573, -0.11016723, -0.11720749,
    1349                          -0.10922083, -0.11748089, -0.10796284,
    1350                          -0.17324244, -0.15860786, -0.16784961,
    1351                          -0.15785693, -0.16501143, -0.1558837,
    1352                          -0.1618709,  -0.15384493, -0.15979744,
    1353                          -0.15251944, -0.15976661, -0.15094149,
    1354                          -0.17055387, -0.1934305, -0.17869999,
    1355                          -0.19402574, -0.17570318, -0.19030006,
    1356                          -0.17114482, -0.18720482, -0.16970896,
    1357                          -0.18836511, -0.17655148, -0.19499397,
    1358                          -0.14139494, -0.14492939, -0.14365011,
    1359                          -0.14969524, -0.14343543, -0.14754936,
    1360                          -0.13944213, -0.14337135, -0.13651514,
    1361                          -0.14183568, -0.13623283, -0.14529803])
     1338                        [-0.02998628, -0.01520652, -0.03043492,
     1339                        -0.0149132, -0.03004706, -0.01476251,
     1340                        -0.0298215, -0.01467976, -0.02988158,
     1341                        -0.01474662, -0.03036161, -0.01442995,
     1342                        -0.07624583, -0.06297061, -0.07733792,
     1343                        -0.06342237, -0.07695439, -0.06289595,
     1344                        -0.07635559, -0.0626065, -0.07633628,
     1345                        -0.06280072, -0.07739632, -0.06386738,
     1346                        -0.12161738, -0.11028239, -0.1223796,
     1347                        -0.11095953, -0.12189744, -0.11048616,
     1348                        -0.12074535, -0.10987605, -0.12014311,
     1349                        -0.10976691, -0.12096859, -0.11087692,
     1350                        -0.16868259, -0.15868061, -0.16801135,
     1351                        -0.1588003, -0.16674343, -0.15813323,
     1352                        -0.16457595, -0.15693826, -0.16281096,
     1353                        -0.15585154, -0.16283873, -0.15540068,
     1354                        -0.17450362, -0.19919913, -0.18062882,
     1355                        -0.19764131, -0.17783111, -0.19407213,
     1356                        -0.1736915, -0.19053624, -0.17228678,
     1357                        -0.19105634, -0.17920133, -0.1968828,
     1358                        -0.14244395, -0.14604641, -0.14473537,
     1359                        -0.1506107, -0.14510055, -0.14919522,
     1360                        -0.14175896, -0.14560798, -0.13911658,
     1361                        -0.14439383, -0.13924047, -0.14829043])
     1362                       
    13621363       
    13631364    def test_bedslope_problem_second_order_one_step(self):
    … …  
    18191820            pass
    18201821
     1822
     1823        #print domain.quantities['level'].centroid_values
     1824
     1825       
    18211826        assert allclose(domain.quantities['level'].centroid_values,
    1822                         [3.96042007e-002, 5.61081828e-002, 4.66578380e-002, 5.73165329e-002,
    1823                          4.72542001e-002, 5.74770060e-002, 4.74459150e-002, 5.77550805e-002,
    1824                          4.80791695e-002, 5.85754074e-002, 4.90681598e-002, 6.02682664e-002,
    1825                          1.16686827e-002, 1.75422685e-002, 1.17014731e-002, 2.15810992e-002,
    1826                          1.30549421e-002, 2.14416681e-002, 1.31212934e-002, 2.15486277e-002,
    1827                          1.34996488e-002, 2.24053139e-002, 1.50195194e-002, 2.22306851e-002,
    1828                          -7.26762170e-003, -1.35071582e-003, -7.88143638e-003,
    1829                          -2.18165245e-003, -7.81749271e-003, -1.06732807e-003,
    1830                          -7.76399231e-003, -1.00580353e-003, -7.81877765e-003,
    1831                          -9.81086203e-004, -7.42500800e-003, -2.41412070e-004, 1.86244453e-001,
    1832                          8.79324341e-002, 1.86232625e-001, 8.78313615e-002, 6.12537452e-002,
    1833                          -3.73125664e-002, -6.37550753e-002, -3.73269705e-002, 6.12145063e-002,
    1834                          8.77700677e-002, 1.86257693e-001, 8.79121535e-002, -4.83328632e-002,
    1835                          1.18336527e-001, -4.83328400e-002, 1.18337005e-001, -4.83328850e-002,
    1836                          -6.65776472e-003, -1.73331646e-001, -1.31654218e-001,
    1837                          -1.73332232e-001, -6.66097985e-003, -4.83323869e-002, 1.18339536e-001,
    1838                          -2.48333331e-001, -2.31666623e-001, -2.48333332e-001,
    1839                          -2.31666628e-001, -2.48333332e-001, -2.31666627e-001,
    1840                          -2.48333330e-001, -2.31666575e-001, -2.48333330e-001,
    1841                          -2.31666597e-001, -2.48333329e-001, -2.31666584e-001,
    1842                          -4.65000000e-001, -3.65000000e-001, -4.65000000e-001,
    1843                          -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
    1844                          -4.65000000e-001, -3.65000000e-001, -4.65000000e-001,
    1845                          -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
    1846                          -5.98333333e-001, -5.81666667e-001, -5.98333333e-001,
    1847                          -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
    1848                          -5.98333333e-001, -5.81666667e-001, -5.98333333e-001,
    1849                          -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
    1850                          -6.48333333e-001, -6.31666667e-001, -6.48333333e-001,
    1851                          -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
    1852                          -6.48333333e-001, -6.31666667e-001, -6.48333333e-001,
    1853                          -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
    1854                          -5.31666667e-001, -5.98333333e-001, -5.31666667e-001,
    1855                          -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
    1856                          -5.31666667e-001, -5.98333333e-001, -5.31666667e-001,
    1857                          -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
    1858                          -4.98333333e-001, -4.81666667e-001, -4.98333333e-001,
    1859                          -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
    1860                          -4.98333333e-001, -4.81666667e-001, -4.98333333e-001,
    1861                          -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
    1862                          -5.48333333e-001, -5.31666667e-001, -5.48333333e-001,
    1863                          -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
    1864                          -5.48333333e-001, -5.31666667e-001, -5.48333333e-001,
    1865                          -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
    1866        
     1827[3.95822638e-002,  5.61022588e-002,  4.66437868e-002,  5.73081011e-002,
     1828 4.72394613e-002,  5.74684939e-002,  4.74309483e-002,  5.77458084e-002,
     1829 4.80628177e-002,  5.85656225e-002,  4.90498542e-002,  6.02609831e-002,
     1830 1.18470315e-002,  1.75136443e-002,  1.18035266e-002,  2.15565695e-002,
     1831 1.31620268e-002,  2.14351640e-002,  1.32351076e-002,  2.15450687e-002,
     1832 1.36414028e-002,  2.24274619e-002,  1.51689511e-002,  2.21789655e-002,
     1833-7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003,
     1834-8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004,
     1835-8.10292716e-003, -3.88584984e-004, -7.35226124e-003,  2.73985295e-004,
     1836 1.86166683e-001,  8.74070369e-002,  1.86166712e-001,  8.74035875e-002,
     1837 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002,
     1838 6.11666725e-002,  8.73846774e-002,  1.86166697e-001,  8.74171550e-002,
     1839-4.83333333e-002,  1.18333333e-001, -4.83333333e-002,  1.18333333e-001,
     1840-4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001,
     1841-1.73333333e-001, -6.66666667e-003, -4.83333333e-002,  1.18333333e-001,
     1842-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
     1843-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
     1844-2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001,
     1845-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
     1846-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
     1847-4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001,
     1848-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
     1849-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
     1850-5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001,
     1851-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
     1852-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
     1853-6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001,
     1854-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
     1855-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
     1856-5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001,
     1857-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
     1858-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
     1859-4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001,
     1860-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
     1861-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001,
     1862-5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001])
     1863
    18671864
    18681865       
Note: See TracChangeset for help on using the changeset viewer.