[ncl-talk] Calculating the area below the curve
Anahita Amiri Farahani
aamir003 at ucr.edu
Tue Dec 17 10:15:54 MST 2019
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/0f2153b3/attachment-0001.html>
More information about the ncl-talk
mailing list