| Re: Help with numerical integration [message #51720] |
Tue, 05 December 2006 01:41 |
Wox
Messages: 184
Registered: August 2006
|
Senior Member |
|
|
On Tue, 05 Dec 2006 10:28:32 +0100, Wox <nomail@hotmail.com> wrote:
> You can for example use IDL's INT_TABULATED where
> X=r and F = ( n(r) * r) / sqrt( n(r)^2 * r^2 - c^2 )
Maybe QROMB or QSIMP are better, I don't know. Check the manual and
try them out :-).
|
|
|
|
| Re: Help with numerical integration [message #51721 is a reply to message #51720] |
Tue, 05 December 2006 01:28  |
Wox
Messages: 184
Registered: August 2006
|
Senior Member |
|
|
On 4 Dec 2006 14:31:47 -0800, "Dave" <Confused.Scientist@gmail.com>
wrote:
> Hi all,
>
> I have never been particularly fond of numerical integration and
> generally do it pretty infrequently these days. Nevertheless, I am
> trying to do a 'quick-and-dirty' atmospheric refraction/ray-trace
> calculation and I'm stumped on the integration. The integral reads:
>
> s = \int_{r1}^{r2} ( n(r) * r * dr ) / sqrt( n(r)^2 - c^2 )
>
> where c is a constant. The trouble I'm having is that n is a function
> of r. Thus, I have a set of discrete points for r:
I'm not sure why "n is a function of r" would be a problem for
numerical integration. As long as you can evaluate it, as you clearly
can, there is no problem. Or am I missing something?
You can for example use IDL's INT_TABULATED where
X=r and F = ( n(r) * r) / sqrt( n(r)^2 * r^2 - c^2 )
|
|
|
|
| Re: Help with numerical integration [message #51727 is a reply to message #51721] |
Mon, 04 December 2006 15:13 |
mmeron
Messages: 44
Registered: October 2003
|
Member |
|
|
In article <1165271507.802034.200410@j72g2000cwa.googlegroups.com>, "Dave" <Confused.Scientist@gmail.com> writes:
> Hi all,
>
> I have never been particularly fond of numerical integration and
> generally do it pretty infrequently these days. Nevertheless, I am
> trying to do a 'quick-and-dirty' atmospheric refraction/ray-trace
> calculation and I'm stumped on the integration. The integral reads:
>
> s = \int_{r1}^{r2} ( n(r) * r * dr ) / sqrt( n(r)^2 - c^2 )
>
> where c is a constant. The trouble I'm having is that n is a function
> of r. Thus, I have a set of discrete points for r:
>
> 0.00000 1.00000 2.00000 3.00000 4.00000
> 5.00000 6.00000
> 7.00000 8.00000 9.00000 10.0000 11.0000
> 12.0000 13.0000
> 14.0000 15.0000 16.0000 17.0000 18.0000
> 19.0000 20.0000
> 21.0000 22.0000 23.0000 24.0000 25.0000
> 26.0000 27.0000
> 28.0000 29.0000 30.0000 31.0000 32.0000
> 33.0000 34.0000
> 35.0000 36.0000 37.0000 38.0000 39.0000
> 40.0000 41.0000
> 42.0000 43.0000 44.0000 45.0000 46.0000
> 47.0000 48.0000
> 49.0000 50.0000 51.0000 52.0000 53.0000
> 54.0000 55.0000
> 56.0000 57.0000 58.0000 59.0000 60.0000
> 61.0000 62.0000
> 63.0000 64.0000 65.0000 66.0000 67.0000
> 68.0000 69.0000
> 70.0000 71.0000 72.0000 73.0000 74.0000
> 75.0000 76.0000
> 77.0000 78.0000 79.0000 80.0000 81.0000
> 82.0000 83.0000
> 84.0000 85.0000 86.0000 87.0000 88.0000
> 89.0000 90.0000
> 91.0000 92.0000 93.0000 94.0000 95.0000
> 96.0000 97.0000
> 98.0000 99.0000 100.000 101.000 102.000
> 103.000 104.000
> 105.000 106.000 107.000 108.000 109.000
> 110.000 111.000
> 112.000 113.000 114.000 115.000 116.000
> 117.000 118.000
> 119.000 120.000
>
> and a set of corresponding points for n(r):
>
> 3210.0997 1915.7633 1031.6773 492.45952
> 253.87699 125.97528
> 60.358510 28.798265 14.395186 7.6860971
> 3.7610915 1.4102413
> 1.0227767 1.0044470 1.0012146 1.0003848
> 1.0001596 1.0001319
> 1.0001193 1.0001084 1.0001147 1.0001217
> 1.0001248 1.0001266
> 1.0001274 1.0001257 1.0001244 1.0001259
> 1.0001277 1.0001302
> 1.0001338 1.0001371 1.0001398 1.0001418
> 1.0001443 1.0001467
> 1.0001496 1.0001525 1.0001557 1.0001591
> 1.0001626 1.0001662
> 1.0001701 1.0001742 1.0001787 1.0001837
> 1.0001883 1.0001915
> 1.0001952 1.0001991 1.0002034 1.0002080
> 1.0002133 1.0002177
> 1.0002231 1.0002281 1.0002309 1.0002338
> 1.0002362 1.0002342
> 1.0002335 1.0002329 1.0002331 1.0002304
> 1.0002273 1.0002229
> 1.0002154 1.0002079 1.0002032 1.0001983
> 1.0001917 1.0001836
> 1.0001748 1.0001659 1.0001570 1.0001472
> 1.0001364 1.0001261
> 1.0001162 1.0001071 1.0000989 1.0000926
> 1.0000869 1.0000819
> 1.0000777 1.0000739 1.0000705 1.0000671
> 1.0000638 1.0000605
> 1.0000570 1.0000535 1.0000501 1.0000466
> 1.0000432 1.0000401
> 1.0000372 1.0000344 1.0000319 1.0000295
> 1.0000272 1.0000251
> 1.0000232 1.0000213 1.0000195 1.0000177
> 1.0000159 1.0000143
> 1.0000127 1.0000112 1.0000097 1.0000083
> 1.0000072 1.0000060
> 1.0000052 1.0000044 1.0000039 1.0000034
> 1.0000030 1.0000026
> 1.0000024
>
> I have three similar integrals to evaluate but I'm pretty lost at the
> moment. Sadly, I suspect this isn't as difficult as I am making it.
>
> Any ideas?
>
Equally spaced data, any variant on Simpson will do well here. I'm
using my own, INTEG, you can find it in Midl_lib on the users
contributions page of RSI. but, as I said, any Simpson will do.
Mati Meron | "When you argue with a fool,
meron@cars.uchicago.edu | chances are he is doing just the same"
|
|
|
|
|
|
comp.lang.idl-pvwave archive
Members
Search
Help
Login
Home
