Dear slytech,
When a formula does not have a closed-form solution, an iterative
approach such as Newton's Method can often be used to arrive at an
arbitrarily precise answer by computing a sequence of successively closer
approximations. Arbitrary precision means that the user can set a margin
of error as small as he likes and then carry out the iteration until it
produces an answer within the margin of error.
In some cases, iterative approximation doesn't work because the formula
converges toward false local solutions or diverges altogether in some
regions of the curve. Fortunately, Allard's Law is proven to converge
toward the true solution for all values that are meaningful in the
physical world. The Aerodrome Meteorological Observing Systems Study Group
(AMOSSG) of the International Civil Aviation Organization describes in
its Manual an efficient way to iteratively solve Allard's Law for the
variable you call R. See page 17 of the following PDF document.
Draft AMOSSG Manual: Appendix B [pdf]
http://www.icao.int/anb/SG/AMOSSG/meetings/amossg4/sn5B.pdf
Note that AMOSSG uses different variable names than yours, but otherwise
gives exactly the same formula for Allard's Law. What you call T, the
visual threshold of illumination, is E_T in the Manual; what you call x,
the extinction coefficient, is sigma (Greek lowercase s) in the Manual;
and R, the visual range or visibility, is called V in the manual.
The name of variable I, for luminous intensity, is unchanged. Although
AMOSSG mentions that I should be 1000 cd when calculating visibility for
aerodrome purposes, you are of course free to use any other value without
affecting the validity of the formula or its iterative solution. Also
note that the variable MOR, which stands for Meteorological Optical Range,
is defined on page 16 as 3/s, so it really doesn't introduce anything new.
I shall quote the relevant passages from page 17 before giving you my code
example. I have made some minor typographical corrections to the text.
Call E_T the visual threshold of illumination and I, the luminous
intensity. With
E_T = I * e^(-s*V) / V^2
and by replacing s by 3/MOR, we get:
V = -MOR/3 * ln(E_T / I * V^2) (1)
[...]
Consider the sequence
V_n = -MOR/3 * ln(E_T / I * V_{n-1}^2) = f(V_{n-1}) .
If this sequence converges, it converges toward V, the visibility
sought.
It can be demonstrated that if V_n is greater than V, then V_{n+1}
will be less than V. The sequence V_n is close to the value of V.
If we take V_0 = MOR and if V_1 is less than V_0, we can conclude
that solution V to equation (1) is less than V_0 = MOR. In this
case, the calculation is not necessary, since the visibility
distance given by the law of Allard is less than the MOR. So the
visibility is equal to the MOR, which is good since the sequence
can diverge in such conditions. However, it is possible to show
that if V_1 > MOR, given that V_0 = MOR, the sequence converges.
The iterative calculation of this sequence can then be performed
until the difference between V_n and V_{n+1} is small in relation
to the value of V_n. For example:
abs(V_n - V_{n-1}) / V_n < 0.01 .
In practice, convergence may be slow. It can be very quickly
accelerated using an intervening variable:
Start with V_0 = MOR and calculate V_1 = f(V_0).
Calculate V_01 = (V_0 + 2*V_1) / 3.
Calculate V_2 = f(V_01) and V_12 = (V_01 + 2*V_2) / 3.
Calculate V_3 = f(V_12) and V_23 = (V_12 + 2*V_3) / 3.
The calculation can be continued so forth, but in practice,
the value of V_23 is very close to the required value of V and
the calculation can be ended at the third iteration.
In the C code below, I use your variable names for T and x, but the
Manual's names for elements in the V_n sequence. We are solving for what
you call R, which the Manual calls V.
double f(double V, double T, double I, double MOR) {
return -MOR/3 * log(T / I * V*V);
}
double compute_visibility(double T, double x, double I) {
double MOR, V_0, V_1, V_01;
if (x == 0) {
printf("error: extinction coefficient must not be zero\n");
return -1.0;
}
V_0 = MOR = 3.0/x;
V_1 = f(V_0, T, I, MOR);
if (V_1 < V_0)
return MOR;
V_01 = (V_0 + 2*V_1) / 3.0;
while (abs(V_1 - V_0) / V_1 >= 0.0001) {
V_0 = V_1;
V_1 = f(V_01, T, I, MOR);
V_01 = (V_01 + 2*V_1) / 3.0;
}
return V_1;
}
The margin of error is currently set to 0.0001, or 0.01%. You may adjust
this value in the test condition of the while loop.
I have enjoyed answering your question. If you have any concerns about
the completeness or accuracy of my research, please advise me through
a Clarification Request and give me a chance to fully meet your needs
before you rate this answer.
Regards,
leapinglizard
Search strategy:
allard law range illumination
://www.google.com/search?hs=gWJ&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+range+illumination&btnG=Search
allard law iterative solution
://www.google.com/search?hs=zqy&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+iterative+solution&btnG=Search
allard law approximation range illumination
://www.google.com/search?hs=fSV&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+approximation+range+illumination&btnG=Search |
Request for Answer Clarification by
slytech-ga
on
24 Apr 2006 02:50 PDT
leapinglizard, thanks for your quick response and sample program. Some questions:
1. I don't quite understand what "convergence" means wrt to the
formula. It would be great if you could provide a few good links to
help me understand the solution better.
2. What is the significance of the "if (V_1 < V_O) return MOR;"
portion of the example? Is it related to the scenario the AMOSS
example gives, i.e. para 4.3.6 mentions "... In this case, the
calculation is not necessary, since the
visibility distance given by the law of allard is less than the MOR".
Or is it related to the solution, in which case the example can only
be used when one is not interested in V when visibility is greater
than MOR.
3. I've tried to run the example provided and compare it with some
other sample data. The results match exactly, however it seems to
break whenever "V_1 < V_O" is true (plus a few other instances). If it
was necessary to calculate V under such conditions, can the example be
improved, or is there an alternative solution?
I've posted the sample data and my results using the example program
you provided at this URL: http://ja.slytech.googlepages.com/home. The
sample data shows the results for V for a number of scenarios, e.g.
when MOR is 10000, 3000, 1000, 300, 100 and 30, for E_t of 10^-4 and
10^-6, and and I of 10000, 1000, 100 and 10.
Thanks!
|
Clarification of Answer by
leapinglizard-ga
on
27 Apr 2006 14:03 PDT
1. I don't quite understand what "convergence" means wrt to the
formula. It would be great if you could provide a few good links to
help me understand the solution better.
Convergence just means that the sequence gets closer and closer to a
single value as you go further along. Some sequences converge for all
values; some converge only for certain values, which is called
conditional convergence. The question of whether a sequence of a given
form converges, and under what conditions, is a thorny one that calls
for deep mathematics. There is, in general, no simple way of deciding
this question.
Mathworld: Convergent series
http://mathworld.wolfram.com/ConvergentSeries.html
Mathworld: Convergence Improvement
http://mathworld.wolfram.com/ConvergenceImprovement.html
Mathworld: Limit
http://mathworld.wolfram.com/Limit.html
2. What is the significance of the "if (V_1 < V_O) return MOR;"
portion of the example? Is it related to the scenario the AMOSS
example gives, i.e. para 4.3.6 mentions "... In this case, the
calculation is not necessary, since the visibility distance given
by the law of Allard is less than the MOR".
Yes, your first hypothesis is correct. My code was a direct
implementation of the AMOSS instructions.
3. I've tried to run the example provided and compare it with some
other sample data. The results match exactly, however it seems to
break whenever "V_1 < V_O" is true (plus a few other instances). If it
was necessary to calculate V under such conditions, can the example be
improved, or is there an alternative solution?
You're quite right about the discrepancies. I may have misunderstood
the AMOSS instructions, or perhaps they are meant for a limited set of
circumstances. I find that if I swap the initial V_0 and V_1 values to
ensure that they are in increasing order, my results agree with those
in the table for all except the six lower-left elements.
4839 / 13400 2653 / 5722 1340 / 2468 572 / 935 247 / 373 94 / 133
-1745 / 8646 1497 / 4090 865 / 1881 409 / 749 188 / 309 75 / 113
-2921 / 4839 -474 / 2653 484 / 1340 265 / 572 134 / 247 57 / 94
-7732 / -1745 -917 / 1497 -174 / 865 150 / 409 86 / 188 41 / 75
My revised code is in the following file.
http://plg.uwaterloo.ca/~mlaszlo/answers/allard.c
I've been thinking of a different way to solve Allard's Law that works
for all solutions within a reasonable range, such as [10^-9, 10^9].
Let me get back to you.
leapinglizard
|
,
0 Comments
)