distance between points on the earth
anal...@hotmail.com
precision for degrees to radians and radius of the earth (by the way,
I have seen values that differ by about 4 miles)
(2) checking argument of acos to be between -1 and 1
Seems to be working OK.
Are the generic sin,cos and acos functions sufficient?
Have functions such as DCOS become unnecessary in f95?
Thanks for any responses.
steve
> (1) Using Double Precision for latitude and longitude and using double
> precision for degrees to radians and radius of the earth (by the way,
> I have seen values that differ by about 4 miles)
>
> (2) checking argument of acos to be between -1 and 1
>
> Seems to be working OK.
>
> Are the generic sin,cos and acos functions sufficient?
Yes. (Of course, you failed to note the accuracy to which you
know your input parameters. If you have more that 17 decimal
digits of precision and you want more than 17 decimal digits,
then DOUBLE PRECISION (assuming IEEE binary64 format) is not
sufficient.)
>
> Have functions such as DCOS become unnecessary in f95?
>
Yes.
glen herrmannsfeldt
> Are the generic sin,cos and acos functions sufficient?
> Have functions such as DCOS become unnecessary in f95?
Functions like DCOS are still part of the standard, but
have few uses. One is that they allow older programs to still
compile without changes. Another is that they can be passed as
actual arguments, but the generics can't. That is rarely useful,
though, as the common usage for passing functions as arguments
of for integration and minimization routines, but integrating
or finding the minima of DCOS isn't very useful.
-- glen
Richard Maine
I'd say no....
But only because it was in f77 that they became superfluous. :-)
They are remnants of f66.
--
Richard Maine | Good judgment comes from experience;
email: last name at domain . net | experience comes from bad judgment.
domain: summertriangle | -- Mark Twain
steve
> steve <kar...@comcast.net> wrote:
> > On Oct 15, 6:08 pm, analys...@hotmail.com wrote:
> > > Have functions such as DCOS become unnecessary in f95?
>
> > Yes.
>
> I'd say no....
>
> But only because it was in f77 that they became superfluous. :-)
>
> They are remnants of f66.
>
For moment, I thought you were disagreeing with me. Thankfully,
the smiley alerted me to your intent. I do note that Glen has point
out that DCOS is needed for passing a specific cosine function to
a subprogram as an actual argument.
--
steve
Richard Maine
> I do note that Glen has point
> out that DCOS is needed for passing a specific cosine function to
> a subprogram as an actual argument.
And I've wanted to do that... let's count... that would be exactly zero
times in over 40 years of programming. Maybe by 50 years I'l hit one,
but I doubt it. The few specific intrinsics you can do that with just
don't turn out to be ones that tend to be the ones that usually want to
get passed as actual arguments in real applications. I suppose I can
imagine passing dcos as an actual argument for something like testing
purposes (say to test an integration routine), but that does't count as
a real application. (If you are numerically integrating dcos in a real
algorithm, we need to talk :-)). And sure, I can invent some artificial
situation, but I've never been in one for real. Glen implied similar
things, though less emphatically than I am doing.
Note that even then, it isn't "needed", but is just a convenience. If
you really wanted the functionality of dcos passed as an actual
argument, you could write your own
function my_dcos(x)
integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
end
which you could then pass as the actual argument. Sure that adds a
little overhead, but considering the overhead already involved in
invoking functions passed as actual arguments, and considering the
number of times you are ever going to want anything like this (zero so
far for me), that's not likely to be much of an issue.
robin
| (1) Using Double Precision for latitude and longitude and using double
| precision for degrees to radians and radius of the earth (by the way,
| I have seen values that differ by about 4 miles)
If you are getting inaccuracies, perghaps you should look for other
causes.
The most likely of those is constants being taken as single precision
where you were assuming that they were double precision.
Note that it is necessary to use the double precision form
of constants such as 1.234D0 (1.234 is taken as single precision),
where double precision is needed.
James Van Buskirk
news:1j7na7d.92wynv1cis7piN%nos...@see.signature...
> Note that even then, it isn't "needed", but is just a convenience. If
> you really wanted the functionality of dcos passed as an actual
> argument, you could write your own
> function my_dcos(x)
> integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
> real(dp), intent(in) :: x
> real(dp) :: my_dcos
> my_dcos = cos(x)
> return
> end
> which you could then pass as the actual argument. Sure that adds a
> little overhead, but considering the overhead already involved in
> invoking functions passed as actual arguments, and considering the
> number of times you are ever going to want anything like this (zero so
> far for me), that's not likely to be much of an issue.
We can force it to be necessary.
C:\gfortran\clf\dcos>type dcos.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
! function my_dcos(x)
elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
subroutine test(fun)
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
real(dp) x(3)
x = [1,2,3]
write(*,*) fun(x)
end subroutine test
end module funcs
program start
use funcs
implicit none
call test(my_dcos)
end program start
C:\gfortran\clf\dcos>gfortran -Wall dcos.f90 -odcos
dcos.f90:33.13:
call test(my_dcos)
1
Error: ELEMENTAL non-INTRINSIC procedure 'my_dcos' is not allowed as an
actual a
rgument at (1)
So we can fix this:
C:\gfortran\clf\dcos>type dcos.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
function my_dcos(x)
! elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
subroutine test(fun)
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
real(dp) x(3)
x = [1,2,3]
write(*,*) fun(x)
end subroutine test
end module funcs
program start
use funcs
implicit none
call test(my_dcos)
end program start
C:\gfortran\clf\dcos>gfortran -Wall dcos.f90 -odcos
C:\gfortran\clf\dcos>dcos
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
But now...
C:\gfortran\clf\dcos>ifort dcos.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
dcos.f90(33) : Error: Procedure argument must be PURE [MY_DCOS]
call test(my_dcos)
-------------^
dcos.f90(33) : Error: The PURE attribute of the associated actual procedure
diff
ers from the PURE attribute of the dummy procedure. [MY_DCOS]
call test(my_dcos)
-------------^
dcos.f90(33) : Error: The ELEMENTAL attribute of the associated actual
procedure
differs from the ELEMENTAL attribute of the dummy procedure. [MY_DCOS]
call test(my_dcos)
-------------^
compilation aborted for dcos.f90 (code 1)
Of course I think this restriction from f95 is silly, but I don't
think it has been changed as of f08, see N1723.pdf, C1234.
Now with DCOS, all this just works:
C:\gfortran\clf\dcos>type dcos.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
function my_dcos(x)
! elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
subroutine test(fun)
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
real(dp) x(3)
x = [1,2,3]
write(*,*) fun(x)
end subroutine test
end module funcs
program start
use funcs
implicit none
intrinsic dcos
! call test(my_dcos)
call test(dcos)
end program start
C:\gfortran\clf\dcos>gfortran -Wall dcos.f90 -odcos
C:\gfortran\clf\dcos>dcos
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
C:\gfortran\clf\dcos>ifort dcos.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
dcos.f90(35) : Error: Procedure argument must be PURE [DCOS]
call test(dcos)
-------------^
compilation aborted for dcos.f90 (code 1)
See? Works perfect, even down to exposing a bug in an admittedly
old version of ifort.
--
write(*,*) transfer((/17.392111325966148d0,6.5794487871554595D-85, &
6.0134700243160014d-154/),(/'x'/)); end
Jan Gerrit Kootstra
I never called a function, only assigned it to a variable.
So is it a bug in ifort or a feature in gfortran?
Kind regards,
Jan Gerrit Kootstra
James Van Buskirk
news:16ac1$4ad8025c$4df909b6$16...@news.chello.nl...
> James,
> I never called a function, only assigned it to a variable.
> So is it a bug in ifort or a feature in gfortran?
Since you have not posted previously in this thread I will assume
that you are talking about my [elided] example, just not perhaps
expressing yourself clearly.
In that case, the example where the DCOS intrinsic was the actual
argument associated with an elemental dummy argument was standard-
conforming. Actually ifort's behavior is a bug in a feature of
ifort. That compiler seems to me to build data structures that
keep better track of whether procedure arguments match their
interfaces. If you don't keep track as well as ifort, you don't
always catch when a dummy argument is PURE and the actual argument
is not. You also don't encounter bugs where the compiler misses
the fact that intrinsic procedures are mostly elemental (with most
of the exceptions being transformational or inquiry) but ifort
accepts that risk because procedure mismatch is one of the
commonest 'hard' bugs to squash and the compiler has been good at
detecting it since the DVF days at least.
Arjan
> (1) Using Double Precision for latitude and longitude and using double
> precision for degrees to radians and radius of the earth (by the way,
> I have seen values that differ by about 4 miles)
> (2) checking argument of acos to be between -1 and 1
1: ACOS is not a good choice for this type of problem.
Try to make a Taylor series around +/- 1. That will fail
since the function is infinitely steep. The same feature
will give numerical imprecision in the neighbourhood of +/- 1.
Solution: Apart from constructing COS(theta) and COS(phi), please
also
construct SIN(theta) and SIN(phi) and then use ATAN2 to get your
angles,
not just ACOS.
2: If you are really interested in geodesics: the earth is not a
sphere.
Get the (free) Proj4 library, compile it and use that to estimate
distances.
That package models the earth as an ellipsoid and knows how to
deal
with distances.
Regards,
Arjan
Tobias Burnus
> elemental function my_dcos(x)
> [...]
> call test(my_dcos)
> 1
> Error: ELEMENTAL non-INTRINSIC procedure 'my_dcos' is not allowed as an
> actual argument at (1)
That error is correct per:
"C1228 (R1221) A nonintrinsic elemental procedure shall not be used as
an actual argument." (F2003)
Same wording is in F2008 (C1233) and F95 (in Section 12.4).
(I find it interesting that NAG f95 allows it as "Extension: ...", I
always thought that it is an extremely strict compiler with almost no
legacy extensions. [ifort allows it a well (it prints a warning with
-stand f03).])
> So we can fix this:
>
> function my_dcos(x)
> [...]
> subroutine test(fun)
> interface
> elemental function fun(x)
>
> dcos.f90(33) : Error: The ELEMENTAL attribute of the associated actual
> procedure differs from the ELEMENTAL attribute of the dummy procedure. [MY_DCOS]
Even without looking at the standard, it makes sense that passing a
non-PURE, non-ELEMENTAL procedure to a PURE and ELEMENTAL dummy is not a
good idea.
The standard has (F2003, in 12.4.1.3 Actual arguments associated with
dummy procedure entities):
"If the interface of the dummy argument is explicit, the characteristics
listed in 12.2 shall be the same for the associated actual argument and
the corresponding dummy argument, except that a pure actual argument may
be associated with a dummy argument that is not pure and an elemental
intrinsic actual procedure may be associated with a dummy procedure"
Thus that "ifort" rejects it is OK. I filled a bug report (PR41724) for
gfortran which does not detect this.
> Now with DCOS, all this just works:
>
> C:\gfortran\clf\dcos>type dcos.f90
> subroutine test(fun)
> interface
> elemental function fun(x)
> [...]
> intrinsic dcos
> call test(dcos)
Is this actually valid or not? If I quote again from "12.4.1.3" but
include the trailing parentheses:
"an elemental intrinsic actual procedure may be associated with a dummy
procedure (which is prohibited from being elemental)."
The "which is prohibited" or as F2008 has it "(which cannot be
elemental).", implies that "interface; elemental function fun" is
already invalid.
Actually, I failed to find this restriction elsewhere in the standard.
As the parenthesis is normative, it should be enough, still - especially
the "cannot" somehow implies that this restriction is listed elsewhere
explicitly.
Tobias
James Van Buskirk
news:4AD8335E...@net-b.de...
I read this as some sort of problem with the standard. At least it
says that it lets you pass DCOS as an actual argument without having
to declare the dummy argument as elemental, which would break codes
dating back to at least f77.
It could be that whoever wrote the quoted paragraph assumed that for
a dummy argument to be elemental was generally prohibited without
actually checking, or that the prohibition existed in a previous
draft but was dropped but not fixed in the quoted paragraph.
I can't see why the prohibition against passing a nonintrinsic
elemental procedure as an actual argument exists anyhow. It seems
to me that elemental procedures are deprecated at birth as a
result, and they seem more useful to me than to deserve such a
fate.
You know, there is a similar restriction regarding procedure
pointers:
C:\gfortran\clf\dcos>type ptr.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
abstract interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
contains
! function my_dcos(x)
elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
end module funcs
program start
use funcs
implicit none
procedure(fun), pointer :: f
real(dp) x(3)
! intrinsic dcos
x = [1,2,3]
f => my_dcos
! f => dcos
write(*,*) f(x)
end program start
C:\gfortran\clf\dcos>gfortran -Wall ptr.f90 -optr
C:\gfortran\clf\dcos>ptr
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
Whoops, gfortran let this one slide!
OK, how about the other versions:
C:\gfortran\clf\dcos>type ptr.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
abstract interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
contains
function my_dcos(x)
! elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
end module funcs
program start
use funcs
implicit none
procedure(fun), pointer :: f
real(dp) x(3)
! intrinsic dcos
x = [1,2,3]
f => my_dcos
! f => dcos
write(*,*) f(x)
end program start
C:\gfortran\clf\dcos>gfortran -Wall ptr.f90 -optr
C:\gfortran\clf\dcos>ptr
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
Now, we agree I think that gfortran should have caught that one as
we shouldn't be able to change non-elemental procedures to syntactically
elemental ones merely by pointing at them. The last version:
C:\gfortran\clf\dcos>type ptr.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
abstract interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
contains
function my_dcos(x)
! elemental function my_dcos(x)
! integer :: dp = kind(1.0d0) !-- Or USE an appropriate module
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
! end
end function
end module funcs
program start
use funcs
implicit none
procedure(fun), pointer :: f
real(dp) x(3)
intrinsic dcos
x = [1,2,3]
! f => my_dcos
f => dcos
write(*,*) f(x)
end program start
C:\gfortran\clf\dcos>gfortran -Wall ptr.f90 -optr
C:\gfortran\clf\dcos>ptr
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
Is there a parenthetical normative clause that appears to forbid this,
too?
Tobias Burnus
> > I read this as some sort of problem with the standard. At least it
> > says that it lets you pass DCOS as an actual argument without having
> > to declare the dummy argument as elemental, which would break codes
> > dating back to at least f77.
It is not completely clear to me why there is/would be a problem.
> > It could be that whoever wrote the quoted paragraph assumed that for
> > a dummy argument to be elemental was generally prohibited without
> > actually checking, or that the prohibition existed in a previous
> > draft but was dropped but not fixed in the quoted paragraph.
Maybe - or I have simply missed to find the right spot.
> > I can't see why the prohibition against passing a nonintrinsic
> > elemental procedure as an actual argument exists anyhow.
Me neither; I think the way elemental procedures are commonly
implemented, it should just work - as ifort and NAG f95 show. (And
gfortran would likely show if the error check were removed.) However,
someone might have had a different implementation in mind
(hypothetically or really existing) were allowing would have lead to
problems.
> > You know, there is a similar restriction regarding procedure
> > pointers:
> >
> > abstract interface
> > elemental function fun(x)
> > [...]
> > elemental function my_dcos(x)
> > [...]
> > procedure(fun), pointer :: f
> > f => my_dcos
> > write(*,*) f(x)
> >
> > Whoops, gfortran let this one slide!
Actually, this program looks OK to me. Do you see anything in the
standard which does not allow this one? (Note, the proc-pointer is not
passed as actual argument but the procedure to which it is associated is
just called.)
> > OK, how about the other versions:
> >
> > abstract interface
> > elemental function fun(x)
> > [...]
> > function my_dcos(x)
> > [...]
> > procedure(fun), pointer :: f
> > f => my_dcos
> >
> > Now, we agree I think that gfortran should have caught that one as
> > we shouldn't be able to change non-elemental procedures to
> > syntactically elemental ones merely by pointing at them.
I agree that it is invalid and that gfortran should have caught it. (As
it is not a constraint, the compiler is not required to check this,
however.)
From the F2003 standard (7.4.2.2 Procedure pointer assignment)
"If proc-pointer-object has an explicit interface, its characteristics
shall be the same as proc-target except that proc-target may be pure
even if proc-pointer-object is not pure and proc-target may be an
elemental intrinsic procedure even if proc-pointer-object is not elemental."
> > The last version:
> >
> > abstract interface
> > elemental function fun(x)
> >
> > procedure(fun), pointer :: f
> > intrinsic dcos
> > f => dcos
Looks OK to me - assuming that "fun" has the same interface as "dcos"
(which seems to be the case). Again, you do not use the proc-pointer as
actual argument thus there should be no problem.
Tobias
Carlie Coats
Single precision arithmetic is normally not adequate for
geo-transformations (i.e.., except on a few platforms like
Cray vector machines).
The Earth is a spheroid, and modeling distances on it requires
elliptic integrals (something that I'd like to avoid, even though
I'm a mathematician myself). Not doing so can lead to errors on the
order of 12KM over distances of 3000 KM, according to an evaluation
I did myself a decade or so ago. (North American continental scale
meteorology plots had easily observable placement errors.)
If you want accurate terrestrial distances, you are better off either
using a package (e.g., "proj4" or "gctp") developed by by serious
geography-programmers, or doing the work "by hand" in a geographic
information system like GRASS.
FWIW -- Carlie Coats
Richard Maine
> On 16 October 2009 James Van Buskirk wrote:
> > > I read this as some sort of problem with the standard. At least it
> > > says that it lets you pass DCOS as an actual argument without having
> > > to declare the dummy argument as elemental, which would break codes
> > > dating back to at least f77.
>
> It is not completely clear to me why there is/would be a problem.
I briefly wondered what this was saying also, but then I concluded that
it was probably just worded a bit hastily. I think he is saying that it
would break codes if it was not allowed to pass dcos as an actual
argument to a nonelemental dummy. That would certainly be true, insomuch
as it was allowed to pass dcos as an actual argument and... well... the
dummy arguments certainly were not elemental as there was no elemental
anything then.
> > > It could be that whoever wrote the quoted paragraph assumed that for
> > > a dummy argument to be elemental was generally prohibited without
> > > actually checking, or that the prohibition existed in a previous
> > > draft but was dropped but not fixed in the quoted paragraph.
>
> Maybe - or I have simply missed to find the right spot.
Yeah. I could have sworn there was such a prohibition, but I can't find
it in a quick skim. I can't say whether my skim missed the prohibition,
that parenthetical mention was intended to be the prohibition (which I'd
sure agree would be a strange way to express it), the prohibition was in
an earlier version of the standard as suggested above, or I am just
plain wrong in thinking there was such a prohibition.
Ah. Don't feel like looking further. It is still bugging me though.
Maybe eventually it will bug me enough that I'll go look a bit more
thoroughly (or maybe I'll get lucky and someone else will find it
first).
> > > I can't see why the prohibition against passing a nonintrinsic
> > > elemental procedure as an actual argument exists anyhow.
>
> Me neither; I think the way elemental procedures are commonly
> implemented, it should just work - as ifort and NAG f95 show. (And
> gfortran would likely show if the error check were removed.) However,
> someone might have had a different implementation in mind
> (hypothetically or really existing) were allowing would have lead to
> problems.
I do recall discussion of multiple possibilities for implementation of
elemental procedures. I don't recall the details, but I do recall that
bits of the standard were drafted so as to allow for multiple
possibilities. In particular, I was pretty sure that allowance for
various implementation strategies went into the restrictions on passing
elementals as arguments... the restrictions that I'm currently having
trouble finding. Seems to me that one set of possible implementation
differences involved whether an elemental procedure looked just like a
non-elemental one, with the calling procedure implicitly putting loops
around the scalar call, or whether the loops were inside of the called
procedure. To my knowledge and recollection (very much second hand, if
not third hand), most or all current compilers put the loops around the
call, but the standard was supposed to accomodate either approach.
James Van Buskirk
> On 16 October 2009 James Van Buskirk wrote:
>> > I read this as some sort of problem with the standard. At least it
>> > says that it lets you pass DCOS as an actual argument without having
>> > to declare the dummy argument as elemental, which would break codes
>> > dating back to at least f77.
> It is not completely clear to me why there is/would be a problem.
Thanks to Richard for decoding the above.
I was thinking in terms of N1601.pdf, C728:
"(R742) The proc-target shall not be a nonintrinsic elemental procedure."
Similar in wording to the actual argument case.
>> > OK, how about the other versions:
>> > abstract interface
>> > elemental function fun(x)
>> > [...]
>> > function my_dcos(x)
>> > [...]
>> > procedure(fun), pointer :: f
>> > f => my_dcos
>> > Now, we agree I think that gfortran should have caught that one as
>> > we shouldn't be able to change non-elemental procedures to
>> > syntactically elemental ones merely by pointing at them.
> I agree that it is invalid and that gfortran should have caught it. (As
> it is not a constraint, the compiler is not required to check this,
> however.)
You can always create situations where the compiler can't compare
interfaces, even if all procedures are module procedures, even in
the same module. That said, if the function invoked had one of
the things that don't make sense for an elemental function, the
results in all likelyhood wouldn't be pretty.
> From the F2003 standard (7.4.2.2 Procedure pointer assignment)
> "If proc-pointer-object has an explicit interface, its characteristics
> shall be the same as proc-target except that proc-target may be pure
> even if proc-pointer-object is not pure and proc-target may be an
> elemental intrinsic procedure even if proc-pointer-object is not
> elemental."
>> > The last version:
>> > abstract interface
>> > elemental function fun(x)
>> > procedure(fun), pointer :: f
>> > intrinsic dcos
>> > f => dcos
> Looks OK to me - assuming that "fun" has the same interface as "dcos"
> (which seems to be the case). Again, you do not use the proc-pointer as
> actual argument thus there should be no problem.
Yes, this I claim to be conforming. Note that in the passage from
7.4.2.2 that you quoted above, there is no question of any hint of
prohibiting this syntax. As a consequence we can construct an example
where a non-intrinsic elemental procedure can be passed to another
procedure. Sorry, no time for examples just now, so you will just
have to use your imagination, as also for the Cray pointer case.
Les Neilson
> precision for degrees to radians and radius of the earth (by the way,
> I have seen values that differ by about 4 miles)
>
Ignoring for the moment the ellipsoid problems.
Radius from center to sea level = R0
Radius from center to summit of Mt Everest = R0 + 5miles approximately
:-)
Les
jfh
> <analys...@hotmail.com> wrote in message
But why ignore the ellipsoid problem? The OP seemed to want better
than 6-figure accuracy.
Assuming a spherical earth gives you 3-figure accuracy but not 4. (Try
calculating the
distance between 2 points on the equator that are 180 deg apart. The
shortest route goes
over a pole: the equatorial route is longer.)
The equatorial radius is 6378140 or 6378137 or 6378136 m according to
3 sources quoted in "Global Earth Physics: A Handbook of Physical
Constants" (T.J.Ahrens, ed., American Geophysical Union 1995), p.8,
and the flattening is 1/298.257 or 1/298.257222. A 4-metre
discrepancy in equatorial radius, not 4 miles. And anyway, shouldn't
high-accuracy geodesics be calculated on the geoid instead of any
ellipsoid?
BTW similar considerations apply to finding the direction from any
point P to Mecca. I do not know whether Muslims use the projection
onto the horizontal plane at P of the straight line through the Earth
to Mecca (ignoring general relativistic amendments to Euclidean
geometry?), or the direction in which the true geodesic on the surface
to Mecca leaves P, or whether they still assume that that geodesic is
the great circle calculated for a spherical earth. The differences are
large if P is close to the antipodes of Mecca, in French Polynesia.
-- John Harper
James Van Buskirk
news:hbb630$882$1...@news.eternal-september.org...
> Yes, this I claim to be conforming. Note that in the passage from
> 7.4.2.2 that you quoted above, there is no question of any hint of
> prohibiting this syntax. As a consequence we can construct an example
> where a non-intrinsic elemental procedure can be passed to another
> procedure. Sorry, no time for examples just now, so you will just
> have to use your imagination, as also for the Cray pointer case.
Testing revealed that I was way overoptimistic about the possibility
of passing a non-intrinsic elemental procedure via a C_FUNPTR --
there are even more barriers here than in the procedure pointer case.
Well, that leaves the Cray pointer test:
C:\gfortran\clf\dcos>type cray.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
! function my_dcos(x)
elemental function my_dcos(x)
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
end function
end module funcs
program start
use funcs
implicit none
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
pointer(f,fun)
real(dp) x(3)
intrinsic dcos
x = [1,2,3]
f = loc(my_dcos)
! f = loc(dcos)
write(*,*) fun(x)
end program start
C:\gfortran\clf\dcos>gfortran -fcray-pointer cray.f90 -ocray
cray.f90:30.11:
f = loc(my_dcos)
1
Error: ELEMENTAL non-INTRINSIC procedure 'my_dcos' is not allowed as an
actual a
rgument at (1)
I don't see why that should be a problem. After all, LOC is an
extension. Maybe that's just a blind spot I have in that I can see
the conflict if LOC was actually going to invoke MY_DCOS, but not if
it's just going to take its address. In that case LOC seems more like
an operator than a function, like operator(=>) (you can't overload
this one though, can you?) but gfortran is quite literal about LOC
being a function, not an operator... See what ifort does:
C:\gfortran\clf\dcos>ifort cray.f90
Intel(R) Fortran Compiler for Intel(R) EM64T-based applications, Version 9.1
Build 20061104
Copyright (C) 1985-2006 Intel Corporation. All rights reserved.
Microsoft (R) Incremental Linker Version 8.00.40310.39
Copyright (C) Microsoft Corporation. All rights reserved.
-out:cray.exe
-subsystem:console
cray.obj
C:\gfortran\clf\dcos>cray
0.540302305868140 -0.416146836547142 -0.989992496600445
The 'B' version, with a non-elemental procedure:
C:\gfortran\clf\dcos>type cray.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
function my_dcos(x)
! elemental function my_dcos(x)
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
end function
end module funcs
program start
use funcs
implicit none
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
pointer(f,fun)
real(dp) x(3)
intrinsic dcos
x = [1,2,3]
f = loc(my_dcos)
! f = loc(dcos)
write(*,*) fun(x)
end program start
C:\gfortran\clf\dcos>gfortran -fcray-pointer cray.f90 -ocray
C:\gfortran\clf\dcos>cray
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
ifort likes this, too. In the general case I don't think the
compiler can tell that there is a mismatch, because we could have
computed the value of f arithmetically. In the intrinsic case:
C:\gfortran\clf\dcos>type cray.f90
module funcs
implicit none
integer, parameter :: dp = kind(1.0d0)
contains
function my_dcos(x)
! elemental function my_dcos(x)
integer, parameter :: dp = kind(1.0d0) !-- Or USE an appropriate module
real(dp), intent(in) :: x
real(dp) :: my_dcos
my_dcos = cos(x)
return
end function
end module funcs
program start
use funcs
implicit none
interface
elemental function fun(x)
import
implicit none
real(dp), intent(in) :: x
real(dp) fun
end function fun
end interface
pointer(f,fun)
real(dp) x(3)
intrinsic dcos
x = [1,2,3]
! f = loc(my_dcos)
f = loc(dcos)
write(*,*) fun(x)
end program start
C:\gfortran\clf\dcos>gfortran -fcray-pointer cray.f90 -ocray
C:\gfortran\clf\dcos>cray
0.54030230586813976 -0.41614683654714241 -0.98999249660044542
works as expected, as does ifort.
anal...@hotmail.com
> <analys...@hotmail.com> wrote in message
Nothing so fancy. Every ZIP code has a latitude and longitude and the
idea is to calculate the distance from one to another.
I am using
FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2)
DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT
ACOSINPUT0 = &
& SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS
(RLAT1/57.2957795130823D0) &
& * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 -
RLONG1/57.2957795130823D0)
ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0)
ZIPDIST1 = 3958.75586574D0 * acos(acosinput)
RETURN
END
All the Zipcodes will be in the US. I would say that one percent
accuracy would be more than sufficient - I just want to avoid any
known pitfalls with violent functions such as acos.
glen herrmannsfeldt
> Nothing so fancy. Every ZIP code has a latitude and longitude and the
> idea is to calculate the distance from one to another.
I don't believe that they are known to 16 significant digits.
> I am using
> FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2)
> DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT
A good use for statement functions!
DOUBLE PRECISION X
SIND(X)=SIN(X/57.2957795130823D0)
COSD(X)=COS(X/57.2957795130823D0)
> ACOSINPUT0 = &
> & SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS
> (RLAT1/57.2957795130823D0) &
> & * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 -
> RLONG1/57.2957795130823D0)
> ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0)
> ZIPDIST1 = 3958.75586574D0 * acos(acosinput)
A few extra significant digits there, too.
Is that at sea level?
(Remember global warming will change sea level.)
> RETURN
> END
> All the Zipcodes will be in the US. I would say that one percent
> accuracy would be more than sufficient - I just want to avoid any
> known pitfalls with violent functions such as acos.
-- glen
Tom Micevski
> (1) Using Double Precision for latitude and longitude and using double
> precision for degrees to radians and radius of the earth (by the way,
> I have seen values that differ by about 4 miles)
here are some values and references for some common spheroidal models:
case ('grs80','gda94') ! Australian ellipsoid
a=6378137.00000_rk
b=6356752.31414_rk
invf=298.257222101_rk
f=0.003352810681225_rk
case ('wgs84','GPS') ! GPS ellipsoid
a=6378137.0_rk
b=6356752.3142_rk
invf=298.257223563_rk
f=one/invf
* http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf (ICSM, Geocentric
Datum of Australia Technical Manual Version 2.3, ISBN 0-9579951-0-5)
* http://earth-info.nga.mil/GandG/publications/tr8350.2/wgs84fin.pdf
(Appendix A) (NIMA Technical Report TR8350.2, Third Edition, Amendment
1, 3 Jan 2000)
>
> (2) checking argument of acos to be between -1 and 1
>
> Seems to be working OK.
>
> Are the generic sin,cos and acos functions sufficient?
>
> Have functions such as DCOS become unnecessary in f95?
>
> Thanks for any responses.
if you're after millimetre accuracy between any 2 points on earth, then
you should use the vincenty method:
Vincenty (1975) Direct and inverse solutions of geodesics on the
ellipsoid with application of nested equations, Survey Review, 23(176),
88-93.
some useful links:
! * http://www.ga.gov.au/geodesy/datums/vincenty_inverse.jsp
! Procedure based on:
! * http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf (Chapter 4)
! * http://www.movable-type.co.uk/scripts/latlong-vincenty.html
! * http://en.wikipedia.org/wiki/Vincenty's_formulae
anal...@hotmail.com
> analys...@hotmail.com wrote:
> > (1) Using Double Precision for latitude and longitude and using double
> > precision for degrees to radians and radius of the earth (by the way,
> > I have seen values that differ by about 4 miles)
>
> here are some values and references for some common spheroidal models:
>
> case ('grs80','gda94') ! Australian ellipsoid
> a=6378137.00000_rk
> b=6356752.31414_rk
> invf=298.257222101_rk
> f=0.003352810681225_rk
> case ('wgs84','GPS') ! GPS ellipsoid
> a=6378137.0_rk
> b=6356752.3142_rk
> invf=298.257223563_rk
> f=one/invf
>
> Datum of Australia Technical Manual Version 2.3, ISBN 0-9579951-0-5)
> (Appendix A) (NIMA Technical Report TR8350.2, Third Edition, Amendment
> 1, 3 Jan 2000)
>
>
>
> > (2) checking argument of acos to be between -1 and 1
>
> > Seems to be working OK.
>
> > Are the generic sin,cos and acos functions sufficient?
>
> > Have functions such as DCOS become unnecessary in f95?
>
> > Thanks for any responses.
>
> if you're after millimetre accuracy between any 2 points on earth, then
> you should use the vincenty method:
>
> Vincenty (1975) Direct and inverse solutions of geodesics on the
> ellipsoid with application of nested equations, Survey Review, 23(176),
> 88-93.
>
> some useful links:
> ! Procedure based on:
> ! *http://www.icsm.gov.au/icsm/gda/gdatm/gdav2.3.pdf(Chapter 4)
> ! *http://www.movable-type.co.uk/scripts/latlong-vincenty.html
> ! *http://en.wikipedia.org/wiki/Vincenty's_formulae
Thanks for the references. Every ZIP code has a latitude and
longitude (I guess of some central point in the area covered by the
ZIP code). I am using the following code:
FUNCTION ZIPDIST1(RLAT1,RLONG1,RLAT2,RLONG2)
DOUBLE PRECISION RLAT1,RLONG1,RLAT2,RLONG2,ACOSINPUT0,ACOSINPUT
ACOSINPUT0 = &
& SIN(RLAT1/57.2957795130823D0) * SIN(RLAT2/57.2957795130823D0) + COS
(RLAT1/57.2957795130823D0) &
& * COS(RLAT2/57.2957795130823D0) * COS(RLONG2/57.2957795130823D0 -
RLONG1/57.2957795130823D0)
ACOSINPUT = MIN(MAX( -1.0D0,ACOSINPUT0),1.0D0)
ZIPDIST1 = 3958.75586574D0 * acos(acosinput)
RETURN
END
All the Zipcodes will be in the US and 1 percent accuracy will be
sufficient. For my purposes, does the Haversine method offer a
significant improvement?
Thanks.
Tom Micevski
looks like your using a great circle/spherical trig distance method
http://www.ga.gov.au/geodesy/datums/distance.jsp
i'm unsure of the accuracy of this method, but the link above gives some
indication (although the points are in australia, not usa). it looks
like your 1% accuracy will likely be easily met.