[ncl-talk] Calculating the area below the curve
Dennis Shea
shea at ucar.edu
Tue Dec 17 12:15:16 MST 2019
See attached
The following was the fix:
if (N%2.ne.0) then ; odd # subintervals
N = N-1 ; <=== added
ans = ans + y(N)*(2*dx(N-1)^2 + 3*dx(N-2)*dx(N-1)) \
/(6*(dx(N-2)+dx(N-1) ) )
ans = ans + y(N-1)*(dx(N-1)^2 + 3*dx(N-1)*dx(N-2)) \
/(6*dx(N-2) )
ans = ans - y(N-2)*dx(N-1)^3 \ ; note - sign
/(6*dx(N-2)*(dx(N-2)+dx(N-1)) )
end if
On Tue, Dec 17, 2019 at 10:16 AM Anahita Amiri Farahani <aamir003 at ucr.edu>
wrote:
> The problem is when there is an odd subinterval.
>
> Thanks,
> Ana
>
> On Tue, Dec 17, 2019 at 10:48 AM Anahita Amiri Farahani <aamir003 at ucr.edu>
> wrote:
>
>> Dear Dennis,
>>
>> I used your function, but I faced some problems when I wanted to
>> calculated for different diameters. Now my data does not have any missing
>> values. here is the code and the data:
>>
>> ;; -*- mode: emacs-lisp; coding: emacs-mule; -*-; -*- mode: ncl;-*-
>>
>> ;***************************************************************************
>>
>> undef("simpson_ne")
>> function simpson_ne(x[*]:numeric, y[*]:numeric)
>> local N, ny, n, dx, ans, ddx
>> begin
>> N = dimsizes(x)
>> ny = dimsizes(y)
>> if (N.ne.ny) then
>> print("simpson_ne: FATAL: x and y are different sizes: nx="+N+"
>> ny="+ny)
>> exit
>> end if
>>
>> dx = x(1:)-x(0:N-2)
>> ans = 0.0
>> do n=1,N-2,2
>> ddx = dx(n)+dx(n-1)
>> ans = ans + y(n) *(dx(n)^3 + dx(n-1)^3 + 3*dx(n)*dx(n-1)*ddx) \
>> ; center (n)
>> /(6*dx(n)*dx(n-1))
>> ans = ans + y(n-1)*(2*dx(n-1)^3 - dx(n)^3 + 3*dx(n)*dx(n-1)^2) \
>> ; left (n-1)
>> /(6*dx(n-1)*ddx)
>> ans = ans + y(n+1)*(2*dx(n)^3 - dx(n-1)^3 + 3*dx(n-1)*dx(n)^2) \
>> ; right (n+1)
>> /(6*dx(n)*ddx )
>> end do
>>
>> if (N%2.ne.0) then ; odd # subintervals
>> ans = ans + y(N)*(2*dx(N-1)^2 + 3*dx(N-2)*dx(N-1)) \
>> /(6*(dx(N-2)+dx(N-1) ) )
>> ans = ans + y(N-1)*(dx(N-1)^2 + 3*dx(N-1)*dx(N-2)) \
>> /(6*dx(N-2) )
>> ans = ans - y(N-2)*dx(N-1)^3 \ ; note - sign
>> /(6*dx(N-2)*(dx(N-2)+dx(N-1)) )
>> end if
>>
>> return(ans)
>> end
>>
>> ;=======================================================
>> ; Simpson Rule for unevenly sampled data
>> ;=======================================================
>> ncol = 4
>> ana = readAsciiTable("test2.csv", ncol, "float", 1)
>>
>> ;---There are no missing values for Freshwater
>>
>> x = ana(:,0) ; Diameter
>> y = ana(:,1) ; Synthetic sea spray
>> z = ana(:,2) ; Fresh water
>> ;N = dimsizes(x) ; 146
>> y at _FillValue = 9999.0 ; seawater has _FillValue
>>
>> simp_frs = simpson_ne(x,z)
>> print("simp_frs="+ sprintf("%15.3f", simp_frs))
>>
>> ;---Use non-missing values of seawater
>>
>> ii = ind(.not.ismissing(y))
>> xi = x(ii)
>> yi = y(ii)
>>
>> simp_sea = simpson_ne(xi,yi)
>> print("simp_sea="+ sprintf("%15.3f", simp_sea))
>>
>> ; 40 to 80 nm including less than 40 nm
>>
>> x1=x(0:40)
>> z1=z(0:40)
>> simp_frs1 = simpson_ne(x1,z1)
>> print("simp_frs40-80 ="+ sprintf("%15.3f", simp_frs1))
>> print(simp_frs1/simp_sea)
>>
>>
>> x2=x(41:60)
>> z2=z(41:60)
>> simp_frs2 = simpson_ne(x2,z2)
>> print("simp_frs80-160 ="+ sprintf("%15.3f", simp_frs2))
>> print(simp_frs2/simp_sea)
>>
>>
>> x3=x(61:79)
>> z3=z(61:79)
>> simp_frs3 = simpson_ne(x3,z3)
>> print("simp_frs160-320 ="+ sprintf("%15.3f", simp_frs3))
>> print(simp_frs3/simp_sea)
>>
>> x4=x(80:97)
>> z4=z(80:97)
>> simp_frs4 = simpson_ne(x4,z4)
>> print("simp_frs320-640 ="+ sprintf("%15.3f", simp_frs4))
>> print(simp_frs4/simp_sea)
>>
>> x5=x(98:107)
>> z5=z(98:107)
>> simp_frs5 = simpson_ne(x5,z5)
>> print("simp_frs640-1280 ="+ sprintf("%15.3f", simp_frs5))
>> print(simp_frs5/simp_sea)
>>
>> x6=x(108:117)
>> z6=z(108:117)
>> simp_frs6 = simpson_ne(x6,z6)
>> print("simp_frs1280-2560 ="+ sprintf("%15.3f", simp_frs6))
>> print(simp_frs6/simp_sea)
>>
>> x7=x(118:126)
>> z7=z(118:126)
>> simp_frs7 = simpson_ne(x7,z7)
>> print("simp_frs2560-5120 ="+ sprintf("%15.3f", simp_frs7))
>> print(simp_frs7/simp_sea)
>>
>> x8=x(127:145)
>> z8=z(127:145)
>> simp_frs8 = simpson_ne(x8,z8)
>> print("simp_frs5120-10240 ="+ sprintf("%15.3f", simp_frs8))
>> print(simp_frs8/simp_sea)
>>
>>
>> The data is:
>> Diameter (nm, x axis) Synthetic Seawater Lake Michigan Freshwater
>> 18.36363636 298.0184898 62.8549
>> 19 379.6318367 91.3129
>> 19.72727273 369.7488571 97.3195
>> 20.45454545 452.7832041 106.884
>> 21.18181818 460.6317551 135.746
>> 21.90909091 490.448102 136.587
>> 22.72727273 577.3782245 135.276
>> 23.54545455 619.9722449 153.029
>> 24.45454545 616.3904286 164.088
>> 25.36363636 675.9878776 157.471
>> 26.27272727 673.3952245 196.833
>> 27.27272727 710.2851224 203.958
>> 28.27272727 757.1509184 192.036
>> 29.27272727 789.4895714 242.472
>> 30.36363636 845.1270816 224.274
>> 31.45454545 865.0478367 236.221
>> 32.63636364 869.1107959 283.161
>> 33.81818182 900.6645918 277.881
>> 35 1001.667918 302.734
>> 36.36363636 990.0728776 320.725
>> 37.63636364 986.8977755 318.618
>> 39 1047.028837 328.856
>> 40.45454545 1024.749082 336.364
>> 41.90909091 1079.294633 381.286
>> 43.45454545 1139.292286 393.798
>> 45.09090909 1174.019571 410.977
>> 46.72727273 1195.079714 416.643
>> 48.45454545 1159.955286 419.43
>> 50.18181818 1231.537 476.518
>> 52.09090909 1283.953653 484.361
>> 54 1333.375143 496.661
>> 55.90909091 1312.178286 492.204
>> 58 1401.140633 534.932
>> 60.09090909 1387.08202 551.7
>> 62.27272727 1408.173286 559.25
>> 64.54545455 1501.003959 557.561
>> 67 1432.459245 530.684
>> 69.45454545 1416.976061 574.684
>> 71.90909091 1517.05602 572.641
>> 74.54545455 1498.982857 550.367
>> 77.36363636 1503.290776 584.655
>> 80.18181818 1525.932449 580.766
>> 83.09090909 1498.583388 585.633
>> 86.09090909 1504.881653 595.263
>> 89.27272727 1589.345755 603.119
>> 92.54545455 1597.619224 639.341
>> 95.90909091 1670.584122 610.778
>> 99.45454545 1647.383571 658.014
>> 103.0909091 1665.38951 660.075
>> 106.9090909 1585.073735 677.112
>> 110.8181818 1555.77951 680.088
>> 114.8181818 1527.127694 713.833
>> 119.0909091 1657.148184 770.754
>> 123.4545455 1697.631816 798.023
>> 127.9090909 1640.725429 804.933
>> 132.6363636 1544.45251 828.813
>> 137.4545455 1541.43 798.361
>> 142.5454545 1583.656667 812.687
>> 147.7272727 1511 784.559
>> 153.1818182 1494.74 819.803
>> 158.8181818 1569.143333 813.677
>> 164.6363636 1494.26 862.247
>> 170.6363636 1468.176667 828.698
>> 176.9090909 1496.57 808.084
>> 183.3636364 1469.826667 778.683
>> 190.0909091 1416.35 768.083
>> 197 1385 810.226
>> 204.2727273 1297.126667 786.788
>> 211.7272727 1256.223333 812.941
>> 219.4545455 1331.04 758.239
>> 227.5454545 1332.116667 744.37
>> 235.9090909 1330.28 704.377
>> 244.5454545 1157.136667 637.229
>> 253.4545455 1213.676667 648.659
>> 262.7272727 1174.203333 648.906
>> 272.3636364 1032.176667 636.226
>> 282.3636364 1012.436667 606.987
>> 292.7272727 1003.206667 616.206
>> 303.4545455 1015.94 611.99
>> 314.5454545 1157.15 582.636
>> 326.0909091 1081.606667 518.618
>> 338 1107.173333 510.237
>> 350.3636364 981.0833333 464.282
>> 363.1818182 957.0533333 443.98
>> 376.5454545 925.8433333 407.32
>> 390.3636364 955.98 461.604
>> 404.6363636 866.6066667 440.411
>> 419.4545455 827.4733333 397.705
>> 434.8181818 748.2333333 382
>> 450.7272727 772.28 318.854
>> 467.2727273 672.1266667 323.109
>> 484.3636364 622.98 292.444
>> 502.0909091 673.5266667 263.288
>> 520.4545455 548.6833333 234.776
>> 539.5454545 503.6366667 218.003
>> 559.3636364 399.383 208.224
>> 579.8181818 384.306 178.02
>> 601.0909091 356.98 167.827
>> 644.3484435 328.19 115.044
>> 692.2197989 301.18 106.942
>> 743.9208627 274.241 92.1581
>> 799.451635 242.377 80.5616
>> 859.7695428 215.253 63.2435
>> 923.917159 193.643 49.151
>> 992.8519107 174.463 36.6713
>> 1066.573798 159.423 26.4571
>> 1146.040248 150.035 20.198
>> 1231.251261 137.426 14.8797
>> 1323.164263 121.248 11.8462
>> 1422.736682 107.02 8.71663
>> 1528.053664 95.7569 8.04464
>> 1641.98749 82.8993 5.81748
>> 1764.53816 69.2672 4.74231
>> 1896.6631 57.6725 4.20472
>> 2038.362312 48.669 3.18714
>> 2190.593223 40.8249 2.84154
>> 2353.355831 32.6918 2.03516
>> 2529.522419 23.9757 1.38237
>> 2718.135559 16.7131 0.748785
>> 2921.110106 11.5676 0.383992
>> 3138.446059 7.72549 0.211196
>> 3373.015701 4.75836 0.211196
>> 3624.81903 2.59956 0.0575989
>> 3894.813474 1.29723 0.0191996
>> 4184.913888 0.619433 0.0383992
>> 4497.992552 0.263469 0
>> 4833.09204 0.100302 0
>> 5194.04206 0.0339347 0
>> 5580.842611 0.010811 0
>> 5997.323403 0.0043044 0
>> 6445.399289 0.000700716 0
>> 6926.027698 0.000100102 0
>> 7443.038336 0.000100102 0
>> 7998.346058 0.000100102 0
>> 8594.823146 0.000200205 0
>> 9236.299309 0 0
>> 9928.519107 0 0
>> 10665.73798 0.000200205 0
>> 11460.40248 0 0
>> 12312.51261 0 0
>> 13231.64263 0 0
>> 14227.36682 0 0
>> 15280.53664 0 0
>> 16419.8749 0 0
>> 17645.3816 0 0
>> 18966.631 0 0
>> Thanks,
>> Ana
>>
>> On Wed, Dec 11, 2019 at 12:32 AM Dennis Shea <shea at ucar.edu> wrote:
>>
>>> I think the numbers derived by Rick are "reasonable".
>>>
>>> I am no sure what should be expected when integrating the Seawater
>>> values. There is a LARGE gap between diameters 187.7 and 710.5. The
>>> seawater values are all missing. To integrate across the seawater values
>>> would mean drawing a straight line between 1244.69 and 384.31
>>>
>>>
>>> Diam Sea Fresh
>>> [SNIP]
>>> 168.5 1241.65 819.803
>>> 174.7 1244.69 813.677
>>> 181.1 *9999 * 862.247
>>> 187.7 *9999 * 828.698
>>> [SNIP]
>>> 710.5 *9999* 135.773
>>> 736.5 *9999 * 107.602
>>> 835 384.306 80.561
>>> 898 356.98 63.243
>>> [SNIP]
>>>
>>> Ricks results:
>>> Fresh: 363163.4
>>> Sea: 1067454
>>>
>>> %> ncl TEST_trap.ncl <== Basically, the trapezoidal
>>> rule [lots of print statements]
>>> (0) trap_fresh= 362866.156
>>> (0) trap_sea = 1125439.000
>>>
>>> %> ncl TEST_simpne.ncl <== Implements simpsons method for
>>> unevenly sampled data in NCL
>>> (0) simp_fresh= 363163.344
>>> (0) simp_sea = 1067453.250
>>>
>>> Again, I am not sure how one would interpret the Seawater result due to
>>> the gap.
>>>
>>> Good Luck
>>>
>>> On Mon, Dec 9, 2019 at 8:54 PM Rick Brownrigg via ncl-talk <
>>> ncl-talk at ucar.edu> wrote:
>>>
>>>> Hi Anahita,
>>>>
>>>> In stepping through the code to where this error message originates, I
>>>> believe it to be due to a bug that was introduced several years ago. Note
>>>> that the examples given on the documentation for simpne() (
>>>> http://ncl.ucar.edu/Document/Functions/Built-in/simpne.shtml) also
>>>> fail with the same message (this function must not get used much). The
>>>> test for "3 non-missings" occurs *after* the call to the underlying fortran
>>>> routine that actually calculates the area. With a debugger, I was able to
>>>> side-step the test, and for your particular dataset, I get:
>>>>
>>>> print(simpne_sea)
>>>> Variable: simpne_sea
>>>> Type: float
>>>> Total Size: 4 bytes
>>>> 1 values
>>>> Number of Dimensions: 1
>>>> Dimensions and sizes: [1]
>>>> Coordinates:
>>>> (0) 1067454
>>>>
>>>> ncl 11> print(simpne_fresh)
>>>> Variable: simpne_fresh
>>>> Type: float
>>>> Total Size: 4 bytes
>>>> 1 values
>>>> Number of Dimensions: 1
>>>> Dimensions and sizes: [1]
>>>> Coordinates:
>>>> (0) 363163.4
>>>>
>>>> Do these values make sense? If I similarly do the "side step" for the
>>>> examples on the docs page, I get the results shown there, so I'm fairly
>>>> confident in the analysis that its the test for non-missings, not the area
>>>> computations, that is erroneous.
>>>>
>>>> Unfortunately, this bug occurs in compiled code. I am no longer in a
>>>> position to i) officially address the bug, and ii) offer you a temporary
>>>> version of NCL until the big-fix appears in the next release. The best I
>>>> can do is file an issue on the github site with as much detail as I know
>>>> now, and maybe a fix will appear in a future release.
>>>>
>>>> I wish I had a better answer.
>>>>
>>>> Rick
>>>>
>>>> On Mon, Dec 9, 2019 at 3:31 PM Anahita Amiri Farahani via ncl-talk <
>>>> ncl-talk at ucar.edu> wrote:
>>>>
>>>>> Hello All,
>>>>>
>>>>> I tried to calculate the area below the curve but I got the error, I
>>>>> used simpne finction. Here is my code and my data:
>>>>>
>>>>> ana = asciiread("test.csv",(/146,3/),"float")
>>>>> ana at _FillValue = 9999
>>>>>
>>>>> xi= ana(:,0) ; Diameter
>>>>> yi= ana(:,1) ;synthetic sea spray
>>>>> zi= ana(:,2) ; Fresh water
>>>>> simpne_sea = simpne(xi,yi)
>>>>> simpne_fresh = simpne(xi,zi)
>>>>>
>>>>> The error that I got is:
>>>>>
>>>>> fatal:simpne: Must have three or more non-missing values.
>>>>>
>>>>> Diameter Seawater Freshwater
>>>>> 20.2 379.07 62.8549
>>>>> 20.9 363.27 91.3129
>>>>> 21.7 455.041 97.3195
>>>>> 22.5 451.004 106.884
>>>>> 23.3 490.321 135.746
>>>>> 24.1 571.821 136.587
>>>>> 25 619.093 135.276
>>>>> 25.9 616.577 153.029
>>>>> 26.9 691.639 164.088
>>>>> 27.9 664.933 157.471
>>>>> 28.9 717.085 196.833
>>>>> 30 769.125 203.958
>>>>> 31.1 784.47 192.036
>>>>> 32.2 871.042 242.472
>>>>> 33.4 881.308 224.274
>>>>> 34.6 871.322 236.221
>>>>> 35.9 902.276 283.161
>>>>> 37.2 1005.14 277.881
>>>>> 38.5 997.639 302.734
>>>>> 40 1007.47 320.725
>>>>> 41.4 1060.78 318.618
>>>>> 42.9 1046.8 328.856
>>>>> 44.5 1095.6 336.364
>>>>> 46.1 1164.61 381.286
>>>>> 47.8 1192.6 393.798
>>>>> 49.6 1226.13 410.977
>>>>> 51.4 1184.22 416.643
>>>>> 53.3 1253.15 419.43
>>>>> 55.2 1297.6 476.518
>>>>> 57.3 1355.38 484.361
>>>>> 59.4 1336.36 496.661
>>>>> 61.5 1409.51 492.204
>>>>> 63.8 1399.77 534.932
>>>>> 66.1 1441.49 551.7
>>>>> 68.5 1537.13 559.25
>>>>> 71 1470.6 557.561
>>>>> 73.7 1454.18 530.684
>>>>> 76.4 1546.17 574.684
>>>>> 79.1 1526.6 572.641
>>>>> 82 1545.55 550.367
>>>>> 85.1 1556.94 584.655
>>>>> 88.2 1527.92 580.766
>>>>> 91.4 1544.07 585.633
>>>>> 94.7 1621.34 595.263
>>>>> 98.2 1628.2 603.119
>>>>> 101.8 1703.88 639.341
>>>>> 105.5 1686.83 610.778
>>>>> 109.4 1679.66 658.014
>>>>> 113.4 1595.34 660.075
>>>>> 117.6 1590.48 677.112
>>>>> 121.9 1554.42 680.088
>>>>> 126.3 1690.79 713.833
>>>>> 131 1714.44 770.754
>>>>> 135.8 1661.12 798.023
>>>>> 140.7 1574.7 804.933
>>>>> 145.9 1463.69 828.813
>>>>> 151.2 1359.23 798.361
>>>>> 156.8 1319.38 812.687
>>>>> 162.5 1231.97 784.559
>>>>> 168.5 1241.65 819.803
>>>>> 174.7 1244.69 813.677
>>>>> 181.1 9999 862.247
>>>>> 187.7 9999 828.698
>>>>> 194.6 9999 808.084
>>>>> 201.7 9999 778.683
>>>>> 209.1 9999 768.083
>>>>> 216.7 9999 810.226
>>>>> 224.7 9999 786.788
>>>>> 232.9 9999 812.941
>>>>> 241.4 9999 758.239
>>>>> 250.3 9999 744.37
>>>>> 259.5 9999 704.377
>>>>> 269 9999 637.229
>>>>> 278.8 9999 648.659
>>>>> 289 9999 648.906
>>>>> 299.6 9999 636.226
>>>>> 310.6 9999 606.987
>>>>> 322 9999 616.206
>>>>> 333.8 9999 611.99
>>>>> 346 9999 582.636
>>>>> 358.7 9999 518.618
>>>>> 371.8 9999 510.237
>>>>> 385.4 9999 464.282
>>>>> 399.5 9999 443.98
>>>>> 414.2 9999 407.32
>>>>> 429.4 9999 461.604
>>>>> 445.1 9999 440.411
>>>>> 461.4 9999 397.705
>>>>> 478.3 9999 382
>>>>> 495.8 9999 318.854
>>>>> 514 9999 323.109
>>>>> 532.8 9999 292.444
>>>>> 552.3 9999 263.288
>>>>> 572.5 9999 234.776
>>>>> 593.5 9999 218.003
>>>>> 615.3 9999 208.224
>>>>> 637.8 9999 178.02
>>>>> 661.2 9999 167.827
>>>>> 685.4 9999 141.841
>>>>> 710.5 9999 135.773
>>>>> 736.5 9999 107.602
>>>>> 835 384.306 80.5616
>>>>> 898 356.98 63.2435
>>>>> 965 328.19 49.151
>>>>> 1037 301.18 36.6713
>>>>> 1114 274.241 26.4571
>>>>> 1197 242.377 20.198
>>>>> 1286 215.253 14.8797
>>>>> 1382 193.643 11.8462
>>>>> 1486 174.463 8.71663
>>>>> 1596 159.423 8.04464
>>>>> 1715 150.035 5.81748
>>>>> 1843 137.426 4.74231
>>>>> 1981 121.248 4.20472
>>>>> 2129 107.02 3.18714
>>>>> 2288 95.7569 2.84154
>>>>> 2458 82.8993 2.03516
>>>>> 2642 69.2672 1.38237
>>>>> 2839 57.6725 0.748785
>>>>> 3051 48.669 0.383992
>>>>> 3278 40.8249 0.211196
>>>>> 3523 32.6918 0.211196
>>>>> 3786 23.9757 0.0575989
>>>>> 4068 16.7131 0.0191996
>>>>> 4371 11.5676 0.0383992
>>>>> 4698 7.72549 0
>>>>> 5048 4.75836 0
>>>>> 5425 2.59956 0
>>>>> 5829 1.29723 0
>>>>> 6264 0.619433 0
>>>>> 6732 0.263469 0
>>>>> 7234 0.100302 0
>>>>> 7774 0.0339347 0
>>>>> 8354 0.010811 0
>>>>> 8977 0.0043044 0
>>>>> 9647 0.00070072 0
>>>>> 10370 0.0001001 0
>>>>> 11140 0.0001001 0
>>>>> 11970 0.0001001 0
>>>>> 12860 0.00020021 0
>>>>> 13820 0 0
>>>>> 14860 0 0
>>>>> 15960 0.00020021 0
>>>>> 17150 0 0
>>>>> 18430 0 0
>>>>> 19810 0 0
>>>>> _______________________________________________
>>>>> ncl-talk mailing list
>>>>> ncl-talk at ucar.edu
>>>>> List instructions, subscriber options, unsubscribe:
>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>
>>>> _______________________________________________
>>>> ncl-talk mailing list
>>>> ncl-talk at ucar.edu
>>>> List instructions, subscriber options, unsubscribe:
>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>
>>>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191217/8ff07009/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: simpson_ne.ncl
Type: application/octet-stream
Size: 1419 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191217/8ff07009/attachment.obj>
More information about the ncl-talk
mailing list