サブルーチンのファイル(distance.f90): ダウンロード
subroutine distance(lat1,lon1,lat2,lon2,d)
!
! Host: aofd30
! Directory: /work2/am/12.Work11/21.Climatology/07.Tend_Adv/src
!
! Revision history:
! 2011-04-05 16:48
! Initial Version
implicit real*8 (a-h,o-z)
real(8),intent(in)::lat1,lon1,lat2,lon2
real(8),intent(out)::d
pi=asin(1.0d0)*2.0d0
! GRS1980
a=6378.137d0
f=1.0d0/298.257222101d0
e2=f*(2.0d0-f)
ed2=e2/(1.0d0-f)**2
phied=lat1*pi/180.0d0
phisd=lat2*pi/180.0d0
alame=lon1*pi/180.0d0
alams=lon2*pi/180.0d0
call geocen(phied,e2,phie)
call geocen(phisd,e2,phis)
cphie = cos(phie)
sphie = sin(phie)
cphis = cos(phis)
sphis = sin(phis)
ae=cphie*cos(alame)
be=cphie*sin(alame)
ce=sphie
as=cphis*cos(alams)
bs=cphis*sin(alams)
cs=sphis
sd2=sqrt((ae-as)**2+(be-bs)**2+(ce-cs)**2)/2.0d0
cd2=sqrt(1.0d0-sd2*sd2)
if(sd2==0.0d0)then
print *
print *,'Error in subroutine distance: sd2=0'
print *,'ae,as,be,bs,ce,cs=',ae,as,be,bs,ce,cs
print *
stop "Terminated."
endif
if(cd2==0.0d0)then
print *
print *,"Error in subroutine distance: cd2=0"
print *
stop "Terminated."
endif
sd=2.0d0*cd2*sd2
diff=alams - alame
20 continue
if(diff.gt.-pi.and.diff.le.pi) then
goto 50
else if(diff.gt.pi) then
diff=diff-2.0d0*pi
else if(diff.le.-pi) then
diff=diff+2.0d0*pi
endif
goto 20
50 continue
phid = (phied+phisd)/2.0d0
sphid = sin(phid)
ro = a * (1.0d0 -e2*sphid**2/2.0d0 &
& +e2*e2*sphid**2/2.0d0 &
& -5.0d0*e2*e2*sphid**4/8.0d0)
d = asin(sd2)*2.0d0 * ro
end subroutine distance
!
subroutine geocen(phid,e2,phi)
implicit real*8 (a-h,o-z)
tphi=(1.0d0-e2)*tan(phid)
phi = atan(tphi)
return
end subroutine geocen
!
! 地球楕円体の長半径と離心率として，GRS1980 の値を用いている。
!
!使用法： call distance (lat1, lon1, lat2, lon2, d)
!
!
!implicit real*8 (a-h,o-z) を宣言しているので，サブルーチン内の実数はすべて倍精度型である．
!変数 説明
!lat1 始点の緯度 (degree) (入力)
!lon1 始点の経度 (degree) (入力)
!lat2 終点の緯度 (degree) (入力)
!lon2 終点の経度 (degree) (入力)
!d 始点から終点までの距離 (km) (出力)
使い方
call distance(lat1,lon1,lat2,lon2,d)
lat1,lon1: 始点の緯度、経度（単位は度）
lat2,lon2: 終点の緯度、経度（単位は度）
d: 距離（単位はkm）
注意
lat1,lon1,lat2,lon2,dはすべて倍精度とする。単精度の変数を使うと間違った結果となる。
参考
単精度→倍精度への変換
b=dble(a)
aを倍精度に変換してbに代入
c=sngl(d)
倍精度の変数dを単精度に変換してcに代入
Perl
http://www.perlmonks.org/?node_id=150054
いろいろ
Calculate distance between two points on a globe
sdist.f90
注意
こちらのサブルーチンに関しては、rlat1,rlon1,rlat2,rlon2はすべて単精度とする。倍精度の変数を使うと間違った結果となる。
! Coded by Hideaki Hase at RIAM, Kyushu Univ. in 1997
! Modified by A. Manda on Jan 15th, 2013
!
! CALCULATE DISTANCE BETWEEN TWO POINTS ON A SPHERE.
!
! REMARK
!
! o(rlon2,rlat2) N
! / ^
! / |
! / S
! o(rlon1,rlat1)
!
!
subroutine sdist(rlat1,rlon1,rlat2,rlon2,dist)
real,intent(in):: rlat1,rlon1,rlat2,rlon2
real,intent(out)::dist
double precision, parameter::PI=3.14159265358979
! rlat1,rlon1 : latitude and longitude of point 1[degree]
! rlat2,rlon2 : latitude and longitude of point 2[degree]
! dist : DISTANCE[km]
double precision x1,x2,y1,y2,coal1,coal2,coal,alpa,stdis &
& ,xrad,yrad,xa,xb,ya,yb,r
double precision, parameter::r=6372.795477598d0 !Earth Radius
xa=dble(rlat1)
xb=dble(rlat2)
ya=dble(rlon1)
yb=dble(rlon2)
! Avoid NaN (Not a Number)
if(xa.eq.xb.and.ya.eq.yb)then
dist = 0.0
return
end if
x1 = xrad(xa)
x2 = xrad(xb)
y1 = yrad(ya)
y2 = yrad(yb)
coal1 = dsin(x1)*dsin(x2)*dcos(y1)*dcos(y2)
coal2 = dsin(x1)*dsin(x2)*dsin(y1)*dsin(y2)+dcos(x1)*dcos(x2)
coal = coal1+coal2
alpa = dacos(coal)
dist = SNGL(r*alpa) ! double -> single
return
end subroutine sdist
double precision function xrad(x)
double precision, intent(in)::x
double precision, parameter::PI=3.14159265358979
xrad=(90.d0-x)/180.d0*PI
return
end function xrad
double precision function yrad(y)
double precision, intent(in)::y
double precision, parameter::PI=3.14159265358979
yrad=y/180.*PI
return
end function yrad
! Earth radius
! http://en.wikipedia.org/wiki/Earth_radius
L=R0*d*(pi/180.0)
L: 距離(km)
d: 緯度の変化量（度）
pi: 円周率(3.141592653589793)
注：この式は経度が異なる場合は使えない。
Fortranの場合には、あらかじめ
real L
と宣言しておかないと正しい結果が得られない。
L=R0*cos(theta*(pi/180.0))*d*(pi/180.0)
L: 距離(km)
R0: 地球の半径≒6378.1 km（赤道半径6378.137 km、極半径が 6356.752 km）
theta: 緯度(度)
d: 経度の変化量（度）
pi: 円周率(3.141592653589793)
注：この式は緯度が異なる場合は使えない。
Fortranの場合には、あらかじめ
real L
と宣言しておかないと正しい結果が得られない。