Re: [PD] haversine formula in Pd

2015-06-08 Thread Cyrille Henry

ahah!!
using -118.4 better than 118.4 also gives me the expected result. (2287.26)

all confusion came that i i was using the wrong set of value...


cheer
c


Le 08/06/2015 11:10, Lorenzo Sutton a écrit :

Hi,

On 07/06/2015 04:48, Max wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Merci Cyrille,

in the formula the intermediate steps are quite small fractions and it
seems their precision is important.
In the test case the Pd implementation is 8917.74 km off the proper
result (2887.26). However I need a precision of about 1m.

So I assume the haversine formula is not implementable in Pd at all?
(unless double precision will be there that is)


I had a go at immplementing it in Pd Vanilla, with a few [expr], and the result 
seems the one expected... no?

Lorenzo


___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list



___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


Re: [PD] haversine formula in Pd

2015-06-08 Thread Lorenzo Sutton

Hi,

On 07/06/2015 04:48, Max wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Merci Cyrille,

in the formula the intermediate steps are quite small fractions and it
seems their precision is important.
In the test case the Pd implementation is 8917.74 km off the proper
result (2887.26). However I need a precision of about 1m.

So I assume the haversine formula is not implementable in Pd at all?
(unless double precision will be there that is)


I had a go at immplementing it in Pd Vanilla, with a few [expr], and the 
result seems the one expected... no?


Lorenzo
#N canvas 314 73 858 610 10;
#X floatatom 168 116 5 0 0 0 - - -, f 5;
#X floatatom 270 115 5 0 0 0 - - -, f 5;
#X floatatom 493 119 8 0 0 0 - - -, f 8;
#X floatatom 576 119 8 0 0 0 - - -, f 8;
#X msg -93 179 6372.8;
#X text 179 91 lat1;
#X obj 168 146 -;
#N canvas 0 50 450 300 dec2rad 0;
#X obj 127 32 inlet dec;
#X obj 155 61 r \$0-pi;
#X obj 127 118 *;
#X obj 155 86 / 180;
#X obj 127 144 outlet;
#X connect 0 0 2 0;
#X connect 1 0 3 0;
#X connect 2 0 4 0;
#X connect 3 0 2 1;
#X restore 168 171 pd dec2rad;
#X obj 670 170 atan;
#X floatatom 670 222 10 0 0 0 - - -, f 10;
#X msg 670 143 1;
#X obj 670 197 * 4;
#X obj 670 115 loadbang;
#X obj 670 245 s \$0-pi;
#X text 284 93 lat2;
#X obj 493 150 -;
#N canvas 0 50 450 300 dec2rad 0;
#X obj 127 32 inlet dec;
#X obj 155 61 r \$0-pi;
#X obj 127 118 *;
#X obj 155 86 / 180;
#X obj 127 144 outlet;
#X connect 0 0 2 0;
#X connect 1 0 3 0;
#X connect 2 0 4 0;
#X connect 3 0 2 1;
#X restore 493 175 pd dec2rad;
#X text 499 95 lon2;
#X text 582 102 lon1;
#X text 227 194 dLat;
#N canvas 0 50 450 300 dec2rad 0;
#X obj 127 32 inlet dec;
#X obj 155 61 r \$0-pi;
#X obj 127 118 *;
#X obj 155 86 / 180;
#X obj 127 144 outlet;
#X connect 0 0 2 0;
#X connect 1 0 3 0;
#X connect 2 0 4 0;
#X connect 3 0 2 1;
#X restore 53 171 pd dec2rad;
#N canvas 0 50 450 300 dec2rad 0;
#X obj 127 32 inlet dec;
#X obj 155 61 r \$0-pi;
#X obj 127 118 *;
#X obj 155 86 / 180;
#X obj 127 144 outlet;
#X connect 0 0 2 0;
#X connect 1 0 3 0;
#X connect 2 0 4 0;
#X connect 3 0 2 1;
#X restore 270 170 pd dec2rad;
#X floatatom 493 201 5 0 0 0 - - -, f 5;
#X text 540 198 dLon;
#X obj 168 267 * 0.5;
#X obj 168 289 sin;
#X obj 168 323 expr pow($f1 \, 2);
#X obj 493 263 * 0.5;
#X obj 493 286 sin;
#X obj 493 323 expr pow($f1 \, 2);
#X obj 286 389 +;
#X obj 318 323 expr cos($f1)*cos($f2);
#X obj 475 352 *;
#X obj 286 452 expr 2*asin(sqrt($f1));
#X obj -92 510 *;
#X floatatom -92 532 10 0 0 0 - - -, f 10;
#X obj 132 92 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1
-1;
#X obj 241 91 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1
-1;
#X obj 463 89 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1
-1;
#X obj 554 95 bng 15 250 50 0 empty empty empty 17 7 0 10 -262144 -1
-1;
#X floatatom 316 197 5 0 0 0 - - -, f 5;
#X msg 168 69 36.12;
#X msg 576 71 -86.67;
#X msg 270 68 33.94;
#X msg 493 71 -118.4;
#X obj -93 25 t b b b b b;
#X text -66 161 R;
#X text 458 452 -- c;
#X floatatom 168 194 8 0 0 0 - - -, f 8;
#X floatatom 38 200 5 0 0 0 - - -, f 5;
#X text -3 154 lat1 = radians(lat1);
#X text 277 149 lat2 = radians(lat2);
#X text 356 422 -- a;
#X floatatom 286 423 8 0 0 0 - - -, f 8;
#X obj -93 86 del 1;
#X text 519 452 c = 2*asin(sqrt(a));
#X obj -93 -19 bng 15 250 50 0 empty empty empty 17 7 0 10 -257985
-1 -1;
#X obj -93 1 t b b;
#X text -72 -22 - test;
#X text 152 305 a = sin(dLat/2)**2 + cos(lat1)*cos(lat2) * sin(dLon/2)**2
, f 81;
#X connect 0 0 6 0;
#X connect 0 0 20 0;
#X connect 1 0 6 1;
#X connect 1 0 21 0;
#X connect 2 0 15 0;
#X connect 3 0 15 1;
#X connect 4 0 34 0;
#X connect 6 0 7 0;
#X connect 7 0 48 0;
#X connect 8 0 11 0;
#X connect 9 0 13 0;
#X connect 10 0 8 0;
#X connect 11 0 9 0;
#X connect 12 0 10 0;
#X connect 15 0 16 0;
#X connect 16 0 22 0;
#X connect 20 0 31 0;
#X connect 20 0 49 0;
#X connect 21 0 31 1;
#X connect 21 0 40 0;
#X connect 22 0 27 0;
#X connect 24 0 25 0;
#X connect 25 0 26 0;
#X connect 26 0 30 0;
#X connect 27 0 28 0;
#X connect 28 0 29 0;
#X connect 29 0 32 1;
#X connect 30 0 53 0;
#X connect 31 0 32 0;
#X connect 32 0 30 1;
#X connect 33 0 34 1;
#X connect 34 0 35 0;
#X connect 36 0 0 0;
#X connect 37 0 1 0;
#X connect 38 0 2 0;
#X connect 39 0 3 0;
#X connect 41 0 0 0;
#X connect 42 0 3 0;
#X connect 43 0 1 0;
#X connect 44 0 2 0;
#X connect 45 0 54 0;
#X connect 45 1 41 0;
#X connect 45 2 43 0;
#X connect 45 3 44 0;
#X connect 45 4 42 0;
#X connect 48 0 24 0;
#X connect 53 0 33 0;
#X connect 54 0 4 0;
#X connect 56 0 57 0;
#X connect 57 0 45 0;
#X connect 57 1 45 0;
___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


Re: [PD] haversine formula in Pd

2015-06-08 Thread Max
Yes, I just realized too. Sorry for the confusion.
So it seems it is possible to implement it in Pd! Even with only 32bit
floats.

Attached all three working implementations and I added also the
equirectangular approximation which should be faster and still accurate
for short distances.

Thank you everyone very much!




On 2015년 06월 08일 18:22, Cyrille Henry wrote:
 ahah!!
 using -118.4 better than 118.4 also gives me the expected result. (2287.26)
 
 all confusion came that i i was using the wrong set of value...
 
 
 cheer
 c
 
 
 Le 08/06/2015 11:10, Lorenzo Sutton a écrit :
 Hi,

 On 07/06/2015 04:48, Max wrote:
 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1

 Merci Cyrille,

 in the formula the intermediate steps are quite small fractions and it
 seems their precision is important.
 In the test case the Pd implementation is 8917.74 km off the proper
 result (2887.26). However I need a precision of about 1m.

 So I assume the haversine formula is not implementable in Pd at all?
 (unless double precision will be there that is)

 I had a go at immplementing it in Pd Vanilla, with a few [expr], and
 the result seems the one expected... no?

 Lorenzo


 ___
 Pd-list@lists.iem.at mailing list
 UNSUBSCRIBE and account-management -
 http://lists.puredata.info/listinfo/pd-list

 
 ___
 Pd-list@lists.iem.at mailing list
 UNSUBSCRIBE and account-management -
 http://lists.puredata.info/listinfo/pd-list

#N canvas 512 304 329 363 10;
#X obj 109 112 inlet deg;
#X obj 109 231 outlet rad;
#X obj 172 177 / 45;
#X obj 172 111 loadbang;
#X obj 109 202 * 0.0174533;
#X text 93 282 rad=deg*pi/180;
#X obj 172 154 atan;
#X msg 172 133 1;
#X text 225 36 part of zexy;
#X text 63 66 convert DEGree to RADiant;
#X connect 0 0 4 0;
#X connect 2 0 4 1;
#X connect 3 0 7 0;
#X connect 4 0 1 0;
#X connect 6 0 2 0;
#X connect 7 0 6 0;
#N canvas 788 1029 1262 734 10;
#X obj 30 25 inlet;
#X obj 28 673 outlet;
#X obj 924 19 inlet;
#X obj 971 19 inlet;
#X obj 978 647 v loca.earth-rad;
#X obj 556 688 outlet;
#X obj 514 78 unpack f f;
#X msg 978 617 bang;
#X floatatom 978 676 15 0 0 0 - - -, f 15;
#X text 836 249 lat1 lon1 lat2 lon2;
#X obj 826 109 unpack f f;
#X obj 826 176 deg2rad;
#X obj 938 140 f \$1;
#X obj 997 144 f \$2;
#X obj 826 87 t l b b;
#X obj 826 636 outlet;
#X obj 826 289 expr acos(sin($f1)*sin($f3)+cos($f1)*cos($f3)*cos($f4-$f2))
;
#X obj 826 367 * 6.3728e+06;
#X text 352 17 3 implementation attempts for the haversine formula
;
#X obj 896 173 deg2rad;
#X obj 951 172 deg2rad;
#X obj 1006 173 deg2rad;
#X text 118 334 http://www.movable-type.co.uk/scripts/latlong.html
;
#X obj 28 427 atan2;
#X obj 28 406 sqrt;
#X obj 62 405 sqrt;
#X obj 28 447 * 2;
#X obj 62 384 -;
#X msg 62 362 1 \$1;
#X obj 28 339 t f f;
#X obj 28 210 sin;
#X obj 83 207 sin;
#X obj 138 206 cos;
#X obj 193 206 cos;
#X obj 28 179 / 2;
#X obj 83 182 / 2;
#X obj 28 74 unpack f f;
#X obj 28 125 - \$1;
#X obj 83 128 - \$2;
#X obj 28 98 t f f;
#X obj 28 301 expr $f1 + ($f2 * $f3 * $f4);
#X obj 28 542 * 6372.8;
#X obj 28 148 deg2rad;
#X obj 83 149 deg2rad;
#X obj 138 148 deg2rad;
#X obj 193 149 deg2rad;
#X obj 29 232 t f f;
#X obj 28 257 *;
#X obj 79 232 t f f;
#X obj 78 257 *;
#X obj 193 120 f \$1;
#X obj 83 98 t f b;
#X obj 641 213 f \$1;
#X obj 505 524 sqrt;
#X obj 556 647 * 6372.8;
#X obj 569 243 deg2rad;
#X obj 498 244 deg2rad;
#X obj 641 245 deg2rad;
#X text 577 266 ph1;
#X text 510 270 th1;
#X obj 568 214 - \$2;
#X obj 696 219 f \$1;
#X obj 696 251 deg2rad;
#X obj 572 103 t f b b, f 8;
#X text 644 273 th2;
#X text 711 276 ph2;
#X obj 502 375 sin;
#X obj 530 373 sin;
#X obj 504 408 -;
#X text 534 410 dz;
#X obj 609 357 cos;
#X obj 582 358 cos;
#X obj 587 383 *;
#X obj 647 360 cos;
#X obj 586 412 -;
#X text 620 412 dx;
#X obj 729 384 sin;
#X obj 696 385 cos;
#X obj 701 415 *;
#X text 734 419 dy;
#X obj 498 310 t f f f;
#X obj 503 431 t f f;
#X obj 502 456 *;
#X obj 585 434 t f f;
#X obj 584 459 *;
#X obj 703 437 t f f;
#X obj 702 462 *;
#X obj 502 486 +;
#X obj 504 508 +;
#X obj 505 544 / 2;
#X obj 505 571 expr asin($f1);
#X obj 510 599 * 2;
#X obj 577 285 t f f;
#X obj 641 294 t f f;
#X connect 0 0 6 0;
#X connect 0 0 14 0;
#X connect 0 0 36 0;
#X connect 2 0 12 1;
#X connect 3 0 13 1;
#X connect 4 0 8 0;
#X connect 6 0 56 0;
#X connect 6 1 63 0;
#X connect 7 0 4 0;
#X connect 10 0 11 0;
#X connect 10 1 19 0;
#X connect 11 0 16 0;
#X connect 12 0 20 0;
#X connect 13 0 21 0;
#X connect 14 0 10 0;
#X connect 14 1 12 0;
#X connect 14 2 13 0;
#X connect 16 0 17 0;
#X connect 17 0 15 0;
#X connect 19 0 16 1;
#X connect 20 0 16 2;
#X connect 21 0 16 3;
#X connect 23 0 26 0;
#X connect 24 0 23 0;
#X connect 25 0 23 1;
#X connect 26 0 41 0;
#X connect 27 0 25 0;
#X connect 28 0 27 0;
#X connect 29 0 24 0;
#X connect 29 1 28 0;
#X connect 30 0 46 0;
#X connect 31 0 48 0;
#X connect 32 0 40 2;
#X connect 33 0 40 3;
#X connect 34 0 30 0;
#X connect 35 0 31 0;
#X connect 36 0 39 0;
#X 

Re: [PD] haversine formula in Pd

2015-06-07 Thread Alexandre Torres Porres
some people were saying expr does internal double precision calculation, is
it true? i dont think so...

2015-06-07 3:10 GMT-03:00 Cyrille Henry c...@chnry.net:



 Le 07/06/2015 04:48, Max a écrit :

 -BEGIN PGP SIGNED MESSAGE-
 Hash: SHA1

 Merci Cyrille,

 in the formula the intermediate steps are quite small fractions and it
 seems their precision is important.
 In the test case the Pd implementation is 8917.74 km off the proper
 result (2887.26). However I need a precision of about 1m.

 i don't know how you get the 2887.26 km proper result : it's not the value
 computed by the website you send (using value in the patch).


 So I assume the haversine formula is not implementable in Pd at all?
 (unless double precision will be there that is)

 Or is there a workaround?


 i'll thing about it.
 cheers

 c



 m.


 On 2015년 06월 07일 05:05, Cyrille Henry wrote:

 hello,

 pow 2 did not like negative number. use t f f and * in order to
 compute the square of a number.

 i correct 2 of your formula according to what i found in your
 link. they both gives the same result : 11805 the same value
 computed with your 2nd link gives 12210. for other value, it also
 look close. so i guess the difference is number precision.


 cheers c



 Le 06/06/2015 20:18, Max a écrit :

 I tried to implement the Haversine formula in Pd which should
 give you the distance in km between two lat long coordinates.

 https://en.wikipedia.org/wiki/Haversine_formula
 http://www.movable-type.co.uk/scripts/latlong.html

 One huge obstacle is the imprecision of the 32bit floats, but
 even without that I can't get the formula to work. I kept 3
 failed implementations in the test-case.

 http://rosettacode.org/wiki/Haversine_formula

 Could someone give me a hand please?



 ___
 Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
 account-management -
 http://lists.puredata.info/listinfo/pd-list



 ___
 Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
 account-management - http://lists.puredata.info/listinfo/pd-list


 -BEGIN PGP SIGNATURE-
 Version: GnuPG v1

 iEYEARECAAYFAlVzsO0ACgkQ3EB7kzgMM6KETQCfRvgOmpqbmyhZaH2Y5hzsRHlC
 XOgAnR+6lqxJfjTj5k7vR747UU8ATrqA
 =pBZv
 -END PGP SIGNATURE-

 ___
 Pd-list@lists.iem.at mailing list
 UNSUBSCRIBE and account-management -
 http://lists.puredata.info/listinfo/pd-list


 ___
 Pd-list@lists.iem.at mailing list
 UNSUBSCRIBE and account-management -
 http://lists.puredata.info/listinfo/pd-list

___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


Re: [PD] haversine formula in Pd

2015-06-06 Thread Max
-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Merci Cyrille,

in the formula the intermediate steps are quite small fractions and it
seems their precision is important.
In the test case the Pd implementation is 8917.74 km off the proper
result (2887.26). However I need a precision of about 1m.

So I assume the haversine formula is not implementable in Pd at all?
(unless double precision will be there that is)

Or is there a workaround?

m.


On 2015년 06월 07일 05:05, Cyrille Henry wrote:
 hello,
 
 pow 2 did not like negative number. use t f f and * in order to
 compute the square of a number.
 
 i correct 2 of your formula according to what i found in your
 link. they both gives the same result : 11805 the same value
 computed with your 2nd link gives 12210. for other value, it also
 look close. so i guess the difference is number precision.
 
 
 cheers c
 
 
 
 Le 06/06/2015 20:18, Max a écrit :
 I tried to implement the Haversine formula in Pd which should
 give you the distance in km between two lat long coordinates.
 
 https://en.wikipedia.org/wiki/Haversine_formula 
 http://www.movable-type.co.uk/scripts/latlong.html
 
 One huge obstacle is the imprecision of the 32bit floats, but
 even without that I can't get the formula to work. I kept 3
 failed implementations in the test-case.
 
 http://rosettacode.org/wiki/Haversine_formula
 
 Could someone give me a hand please?
 
 
 
 ___ 
 Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
 account-management - 
 http://lists.puredata.info/listinfo/pd-list
 
 
 
 ___ 
 Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
 account-management - http://lists.puredata.info/listinfo/pd-list
 

-BEGIN PGP SIGNATURE-
Version: GnuPG v1

iEYEARECAAYFAlVzsO0ACgkQ3EB7kzgMM6KETQCfRvgOmpqbmyhZaH2Y5hzsRHlC
XOgAnR+6lqxJfjTj5k7vR747UU8ATrqA
=pBZv
-END PGP SIGNATURE-

___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


Re: [PD] haversine formula in Pd

2015-06-06 Thread Claude Heiland-Allen

On 07/06/15 03:48, Max wrote:

In the test case the Pd implementation is 8917.74 km off the proper
result (2887.26). However I need a precision of about 1m.


Circumference of the Earth in meters: 40,075,000
Accuracy of single precision (24bit): 16,777,216

So your input values will already be inaccurate before you even start 
messing with rounding errors etc.



So I assume the haversine formula is not implementable in Pd at all?
(unless double precision will be there that is)


I guess so.


Or is there a workaround?


Not easy at all, but:

Maybe use two floats at different scales 'a' (close to the true value), 
'b' (the residual difference from the true value) to express each 
coordinate 'c', where:


c = a + b

This would give around 48bits, which should be enough.  Actually 
performing the addition would give just 'a', due to limited precision. 
You have to work with the unpleasant properties of floating point 
numbers, like (a + b) - b not always being equal to a.


There are a few blog posts out there about using it on GPU with GLSL 
(emulated double, double-single), and there is the QD package for 
double-double and quad-double in C++ and FORTRAN (maybe with C 
wrappers), there's a Haskell library called compensated which might 
give some ideas also.


You'll probably also need to do some algebraic manipulations with 
trigonometric identities etc to get accurate results.


See:
https://www.thasler.com/blog/blog/glsl-part2-emu
http://crd-legacy.lbl.gov/~dhbailey/mpdist/
http://hackage.haskell.org/package/compensated
https://en.wikipedia.org/wiki/List_of_trigonometric_identities


Claude
--
http://mathr.co.uk


___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


Re: [PD] haversine formula in Pd

2015-06-06 Thread Joel Matthys
You'll need to do the calculations in an external, where you can use 
double precision.


And. TADA! Here's the external. I just copied in the C code from the 
Rosettacode link you sent. It gives 2887.26 as the result now.


The external code is here: https://github.com/jwmatthys/haversine-pd

I built the linux version already. It should be easy to build others.

Joel

On 06/06/2015 09:48 PM, Max wrote:

-BEGIN PGP SIGNED MESSAGE-
Hash: SHA1

Merci Cyrille,

in the formula the intermediate steps are quite small fractions and it
seems their precision is important.
In the test case the Pd implementation is 8917.74 km off the proper
result (2887.26). However I need a precision of about 1m.

So I assume the haversine formula is not implementable in Pd at all?
(unless double precision will be there that is)

Or is there a workaround?

m.


On 2015년 06월 07일 05:05, Cyrille Henry wrote:

hello,

pow 2 did not like negative number. use t f f and * in order to
compute the square of a number.

i correct 2 of your formula according to what i found in your
link. they both gives the same result : 11805 the same value
computed with your 2nd link gives 12210. for other value, it also
look close. so i guess the difference is number precision.


cheers c



Le 06/06/2015 20:18, Max a écrit :

I tried to implement the Haversine formula in Pd which should
give you the distance in km between two lat long coordinates.

https://en.wikipedia.org/wiki/Haversine_formula
http://www.movable-type.co.uk/scripts/latlong.html

One huge obstacle is the imprecision of the 32bit floats, but
even without that I can't get the formula to work. I kept 3
failed implementations in the test-case.

http://rosettacode.org/wiki/Haversine_formula

Could someone give me a hand please?



___
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
account-management -
http://lists.puredata.info/listinfo/pd-list



___
Pd-list@lists.iem.at mailing list UNSUBSCRIBE and
account-management - http://lists.puredata.info/listinfo/pd-list


-BEGIN PGP SIGNATURE-
Version: GnuPG v1

iEYEARECAAYFAlVzsO0ACgkQ3EB7kzgMM6KETQCfRvgOmpqbmyhZaH2Y5hzsRHlC
XOgAnR+6lqxJfjTj5k7vR747UU8ATrqA
=pBZv
-END PGP SIGNATURE-

___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list



___
Pd-list@lists.iem.at mailing list
UNSUBSCRIBE and account-management - 
http://lists.puredata.info/listinfo/pd-list


[PD] haversine formula in Pd

2015-06-06 Thread Max
I tried to implement the Haversine formula in Pd which should give you
the distance in km between two lat long coordinates.

https://en.wikipedia.org/wiki/Haversine_formula
http://www.movable-type.co.uk/scripts/latlong.html

One huge obstacle is the imprecision of the 32bit floats, but even
without that I can't get the formula to work. I kept 3 failed
implementations in the test-case.

http://rosettacode.org/wiki/Haversine_formula

Could someone give me a hand please?
#N canvas 512 301 329 363 10;
#X obj 109 112 inlet deg;
#X obj 109 231 outlet rad;
#X obj 172 177 / 45;
#X obj 172 111 loadbang;
#X obj 109 202 * 0.0174533;
#X text 93 282 rad=deg*pi/180;
#X obj 172 154 atan;
#X msg 172 133 1;
#X text 225 36 part of zexy;
#X text 63 66 convert DEGree to RADiant;
#X connect 0 0 4 0;
#X connect 2 0 4 1;
#X connect 3 0 7 0;
#X connect 4 0 1 0;
#X connect 6 0 2 0;
#X connect 7 0 6 0;
#N canvas 543 85 1262 734 10;
#X obj 30 25 inlet;
#X text 118 294 http://www.movable-type.co.uk/scripts/latlong.html
;
#X obj 28 387 atan2;
#X obj 28 366 sqrt;
#X obj 62 365 sqrt;
#X obj 28 407 * 2;
#X obj 62 344 -;
#X msg 62 322 1 \$1;
#X obj 28 299 t f f;
#X obj 28 673 outlet;
#X obj 28 210 sin;
#X obj 83 207 sin;
#X obj 138 206 cos;
#X obj 193 206 cos;
#X obj 28 179 / 2;
#X obj 83 182 / 2;
#X obj 28 74 unpack f f;
#X obj 28 125 - \$1;
#X obj 83 128 - \$2;
#X obj 83 98 t f f;
#X obj 28 98 t f f;
#X obj 28 234 pow 2;
#X obj 83 235 pow 2;
#X obj 28 261 expr $f1 + ($f2 * $f3 * $f4);
#X obj 924 19 inlet;
#X obj 971 19 inlet;
#X obj 978 647 v loca.earth-rad;
#X obj 563 625 outlet;
#X obj 514 78 unpack f f;
#X obj 512 185 -;
#X obj 513 133 f \$1;
#X obj 512 157 pack f f;
#X obj 571 189 -;
#X obj 571 161 pack f f;
#X obj 572 137 f \$2;
#X obj 514 105 t b f f;
#X obj 573 109 t b f b;
#X obj 670 161 f \$1;
#X obj 510 251 / 2;
#X obj 510 276 sin;
#X obj 626 257 cos;
#X obj 677 255 cos;
#X obj 568 251 / 2;
#X obj 570 273 sin;
#X obj 632 291 *;
#X obj 595 319 *;
#X obj 595 364 pow 2;
#X obj 562 499 expr asin($f1);
#X obj 561 412 +;
#X obj 562 471 sqrt;
#X obj 562 525 * 2;
#X msg 978 617 bang;
#X floatatom 978 676 15 0 0 0 - - -, f 15;
#X obj 464 323 expr pow($f1 \, 2);
#X obj 469 355 * -1;
#X obj 563 596 * 6372.8;
#X text 836 249 lat1 lon1 lat2 lon2;
#X obj 826 109 unpack f f;
#X obj 826 176 deg2rad;
#X obj 938 140 f \$1;
#X obj 997 144 f \$2;
#X obj 826 87 t l b b;
#X obj 826 636 outlet;
#X obj 826 289 expr acos(sin($f1)*sin($f3)+cos($f1)*cos($f3)*cos($f4-$f2))
;
#X obj 28 502 * 6372.8;
#X obj 826 367 * 6.3728e+06;
#X text 352 17 3 implementation attempts for the haversine formula
;
#X obj 28 148 deg2rad;
#X obj 83 149 deg2rad;
#X obj 138 148 deg2rad;
#X obj 193 149 deg2rad;
#X obj 513 214 deg2rad;
#X obj 568 215 deg2rad;
#X obj 623 214 deg2rad;
#X obj 678 215 deg2rad;
#X obj 896 173 deg2rad;
#X obj 951 172 deg2rad;
#X obj 1006 173 deg2rad;
#X connect 0 0 16 0;
#X connect 0 0 28 0;
#X connect 0 0 61 0;
#X connect 2 0 5 0;
#X connect 3 0 2 0;
#X connect 4 0 2 1;
#X connect 5 0 64 0;
#X connect 6 0 4 0;
#X connect 7 0 6 0;
#X connect 8 0 3 0;
#X connect 8 1 7 0;
#X connect 10 0 21 0;
#X connect 11 0 22 0;
#X connect 12 0 23 2;
#X connect 13 0 23 3;
#X connect 14 0 10 0;
#X connect 15 0 11 0;
#X connect 16 0 20 0;
#X connect 16 1 19 0;
#X connect 17 0 67 0;
#X connect 18 0 68 0;
#X connect 19 0 18 0;
#X connect 19 1 70 0;
#X connect 20 0 17 0;
#X connect 20 1 69 0;
#X connect 21 0 23 0;
#X connect 22 0 23 1;
#X connect 23 0 8 0;
#X connect 24 0 17 1;
#X connect 24 0 30 1;
#X connect 24 0 59 1;
#X connect 25 0 18 1;
#X connect 25 0 34 1;
#X connect 25 0 60 1;
#X connect 26 0 52 0;
#X connect 28 0 35 0;
#X connect 28 1 36 0;
#X connect 29 0 71 0;
#X connect 30 0 31 0;
#X connect 31 0 29 0;
#X connect 32 0 72 0;
#X connect 33 0 32 0;
#X connect 34 0 33 0;
#X connect 35 0 30 0;
#X connect 35 1 31 1;
#X connect 35 2 73 0;
#X connect 36 0 34 0;
#X connect 36 1 33 1;
#X connect 36 2 37 0;
#X connect 37 0 74 0;
#X connect 38 0 39 0;
#X connect 39 0 53 0;
#X connect 40 0 44 0;
#X connect 41 0 44 1;
#X connect 42 0 43 0;
#X connect 43 0 45 0;
#X connect 44 0 45 1;
#X connect 45 0 46 0;
#X connect 46 0 48 1;
#X connect 47 0 50 0;
#X connect 48 0 49 0;
#X connect 49 0 47 0;
#X connect 50 0 55 0;
#X connect 51 0 26 0;
#X connect 53 0 54 0;
#X connect 54 0 48 0;
#X connect 55 0 27 0;
#X connect 57 0 58 0;
#X connect 57 1 75 0;
#X connect 58 0 63 0;
#X connect 59 0 76 0;
#X connect 60 0 77 0;
#X connect 61 0 57 0;
#X connect 61 1 59 0;
#X connect 61 2 60 0;
#X connect 63 0 65 0;
#X connect 64 0 9 0;
#X connect 65 0 62 0;
#X connect 67 0 14 0;
#X connect 68 0 15 0;
#X connect 69 0 12 0;
#X connect 70 0 13 0;
#X connect 71 0 38 0;
#X connect 72 0 42 0;
#X connect 73 0 40 0;
#X connect 74 0 41 0;
#X connect 75 0 63 1;
#X connect 76 0 63 2;
#X connect 77 0 63 3;
#N canvas 50 388 487 549 10;
#X text 248 209 calculates distance of two points;
#X text 66 383 http://rosettacode.org/wiki/Haversine_formula;
#X text 66 328 In this test-case \, let's calculate