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