4 replies to this topic

[Haldis]

Hi, I am trying to make a program on my calc (casio fx 9750g plus) that converts coordinates from ECEF(Earth-Centered, Earth-Fixed) to geodetic coordinates(latitude, longitude) but I just cant get it to work. In the formula i need to iterate, and that is were the problem start. The answer will go towards the solution after a number of loops, but I dont know how to make the formula for a loop on my calc.

This is how the command look in MATLAB:

function [lat,lon,h]=ECEF2geod(a,b,X,Y,Z)

lon=atan(Y/X);
lat_old=0;
ro=sqrt(X^2+Y^2);

for i=1:5
lat_new=atan(Z/ro+((NRad(a,b,lat_old)*e2(a,)*sin(lat_old)/ro));
lat_old=lat_new ;
end
lat=lat_new;

h=ro*cos(lat_new)+Z*sin(lat_new)-NRad(a,b,lat_new)*(1-e2(a,*(sin(lat_new)^2));

Here the loop is set to go 5 times, and it will give a pretty accurate answer. But i do not know how i can make the same program for my calc. Anyone got any tips on how i should do this?

Thanks in advance, and sorry for my bad english.

[Somelauw]

Well, I would recommend to start learning the language.
Look in the manual or search for a good online tutorial.

[verena]

I dont know how to make the formula for a loop on my calc.

Your english is not bad at all, it's just that a loop does not require a "formula" but "program control commands". As somelauw said, look in your manual, at "programing".

There are several ways to program a loop. The most "oldfashioned" one first:

1) set a label for instance "LBL 0" (means: label number zero)
2) calculate the new values for h, lat, long
3) test if the new solution is good enough, or still off in which case you need to repeat the iteration.
The testing is usually done by remembering one or more of the variables calculated, and comparing them to the solutions(s) of the repeat calculation, by using the absolute value of the difference, or its square.

Example: You could for instance test if the "new" lat "lat_new" is 1/100 arc second or less "off" from the "previous" lat ("lat_old"). Use variable J and K to remember and calculate things. And let's call "lat_new" "L" for simplicity, and "lat_old" K.
Assuming you calculate in "degrees" (not "radians")

[Begin of program]
Assume "lat_old" is zero, to begin with. Store 0(zero) in K. :0->K:
(where "->" is the "store" symbol)
[begin of loop]
LBL 0:
(Calculate new lat (L) and long)
Then, calculate abs(L-K) and store in J. Then store L in K, to remember it as "previous result", like in your example: lat_old := lat_new, in casiospeak: :L->K:
Then test:
J>1/360000=>GOTO 0:
(means: if the new L is more than 1/100 arc second off the old h, repeat the calulation from LBL 0)
[end of loop]
NOTE1 ">" is called "relational operator" see your manual, to compare J and 1/360000.
NOTE2 "=>" is the "(conditional) jump" command symbol, consult your manual in the chapter about "programming".
(Display calculated results.)
[end of program]

Well of course the first time through the loop, the difference between zero and L is large, so the program will go back to LBL 0, and recalculate the coordinates.
As the new longitudes get calculated, they approach more and more the value remembered in K. When the difference is 1/360000 degree or less, the program will not go to LBL 0 anymore and the iterations are done.

You can also program loops with "DO" and "LPWHILE":

(store zero in K)
DO:
(calculations for lat, long, J, and store L in K)
LPWHILE J>1/360000:
Here, DO is like the lable "LBL 0" and LPWHILE functions like "IF J larger than 1/360000, GOTO DO"

or "WHILE" and "WHILEEND", are all slight variations to the same theme:

(store 1 in J)
WHILE J>1/360000:
(calculations etc.)
WHILEEND:
Here, the test is at the beginning, WHILEEND makes program control go automatically back to WHILE, and if the test at the WHILE fails (when J is equal or smaller than 1/360000), the program goes to the first statement after the "WHILEEND".

Have fun

[Haldis]

Thanks a lot. Now i understand a little more how the loop function works. I still cant get it to give me the answer i want. The longitude is ok, but the latitude is way of. Here is how the program look in my calc (only for latitude as the higth and longitude is easy):

"LATOLD"?->K:"A"?->A:"B"?->B:"X"?->X:"Y"?->Y:"Z")->Z:
(X^2+Y^2)->R
((A^2-B^2)/A^2)->E
Do:A/sqrt(1-(E*(sin K)^2))->N
"LAT"_ (_ is equal to the output triangle)
atan(Z/R+(N*E(sin K))/R)->L
Abs(L-K)->J
L->K
LpWhile J>1/360000:

How do it look? It migth be a error in my formula but i think it is because it does not loop properly.

[verena]

Well I only tried to explain how a loop works...

1) You could try to make the loop go five times for instance by introducing a counter:

before the loop, program: 5 -> C

in the loop, program C - 1 -> C

and test with LPWHILE C > 0

Then see if it's still off after five loops

2) What is the purpose of the a and b in the geodetic formula ? And what do NRad and E2 do ? (Also, it doesn't help that your original formula now has smileys in it, please fix)

3) Did you set the calc to calculate angles in degrees ? It's good practice to put that inside your prgram as well.