RTKLIB v2.4.2が公開されていた。
毎日HPを確認していたんだが、ついに！
http://www.rtklib.com/rtklib.htm
ベータ版はほとんど見てないけど、v2.4.1からどれくらい変わったのだろうか？
PPP-ARの仕組みは全く勉強していないので、これを機にソースを読んで勉強したい。
2013年4月30日火曜日
2013年4月29日月曜日
Okada(1992) DC3Dの使い方
Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040]の式のfortran77プログラムが下記で公開されている。
http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html
ある場所（緯度経度）での変位を出力しようと思ったら、結構頭がこんがらがった。というわけで、wrapperを作ったので、プログラムをメモ代わりにアップしておく。断層パラメータの記載は一般的な定義のままで、位置は断層上端左端にした。ああ、ややこしかった。
注意としては、深さは正の値を入れること、長さLの場合はAL1=0, AL2=Lとすること、幅Wの場合は、AW1=-W, AW2=0とすること。
一般的な断層パラメータの定義が数式と少し違い（走向は北から反時計回りとか）、注意が必要ですね。すぐ忘れそうなので、やはりwrapperを作っておいた方が良いな。もう少しモジュールがたまってきたら、arでライブラリにしてしまいたい。
ちなみに、こんな感じで呼び出す。
モジュールは下記のようなもの。
http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_J.html
ある場所（緯度経度）での変位を出力しようと思ったら、結構頭がこんがらがった。というわけで、wrapperを作ったので、プログラムをメモ代わりにアップしておく。断層パラメータの記載は一般的な定義のままで、位置は断層上端左端にした。ああ、ややこしかった。
注意としては、深さは正の値を入れること、長さLの場合はAL1=0, AL2=Lとすること、幅Wの場合は、AW1=-W, AW2=0とすること。
一般的な断層パラメータの定義が数式と少し違い（走向は北から反時計回りとか）、注意が必要ですね。すぐ忘れそうなので、やはりwrapperを作っておいた方が良いな。もう少しモジュールがたまってきたら、arでライブラリにしてしまいたい。
ちなみに、こんな感じで呼び出す。
test.f90:
program test
use okada92_util
implicit none
type(disp) :: du
type(flt) :: f
real(8) :: lat, lon, h
integer :: i, j
f%lambda = 1.d0 ! ずるしている
f%mu = 1.d0
f%lat = 32.7725d0
f%lon = 136.0282d0
f%depth = 5000.d0
f%str = 261.766d0
f%dip = 8.903d0
f%rake = 90.d0
!f%rake = 0.d0
f%length = 100000.d0
f%width = 25263.533d0
f%slip = 5.d0
do i = -50, 50
do j = -50, 50
lat = f%lat + dble(i)*0.1d0
lon = f%lon + dble(j)*0.1d0
h = 0.d0
du = flt_disp(f, lat, lon, h)
write(*,'(F7.3,",",F7.3,3(",",F15.7))') lon, lat, du%e, du%n, du%u
enddo
enddo
stop
end program test
モジュールは下記のようなもの。
okada92_util.f90:
module okada92_util
type disp
real(8) :: e, n, u
end type disp
type flt
!real(8) :: alpha ! (lambda+mu)/(lambda+2*mu) !αをインプットにした方が良かったかも
real(8) :: lambda, mu ! (lambda+mu)/(lambda+2*mu)
real(8) :: lat, lon ! latitude, longitude (deg) at the corner of upper edge
real(8) :: depth ! fault depth (positive is ...)
real(8) :: str ! strike angle (deg)
real(8) :: dip ! dip angle (deg)
real(8) :: rake ! rake angle (deg)
real(8) :: length, width ! fault length, width (m)
real(8) :: slip ! slip (m)
end type flt
! ellipsoid
type ellip
real(8) :: a, b, e2, f, fi
end type ellip
! parameters
real(8), parameter :: PI = 3.1415926535897932d0
real(8), parameter :: D2R = PI/180.d0
real(8), parameter :: R2D = 180.d0/PI
contains
! ----------------------------------------------------------
type(disp) function flt_disp(f, lat, lon, h)
! ----------------------------------------------------------
! Calculates displacement at (lat, lon, h) due to a
! rectangular fault. The equations are given by
! Okada (1992) [Bull. Seism. Soc. Am., 82, 1018-1040].
!
! This function is a wrapper of 'DC3D' on the website:
! http://www.bosai.go.jp/study/application/dc3d/DC3Dhtml_E.html
!
! Input:
! f : type(flt), includes
! f%lambda, mu ! (lambda+mu)/(lambda+2*mu)
! f%lat, lon ! latitude, longitude (deg) at the corner of upper edge.
! f%depth ! fault depth (positive is downward)
! f%str ! strike angle (deg)
! f%dip ! dip angle (deg)
! f%rake ! rake angle (deg)
! f%length, width ! fault length, width (m)
! f%slip ! slip (m)
! lat, lon, h : latitude(deg), longitude(deg), h(m)
! positive h is upward.
! Output:
! disp : type(disp),
! disp%e ! eastward disp. (m)
! disp%n ! northward disp. (m)
! disp%u ! vertical disp. (m)
!
! Created:
! 2013/04/29 Sako
! ----------------------------------------------------------
type(flt), intent(in) :: f
real(8), intent(in) :: lat, lon, h
! inputs of DC3D
real(8) :: a, x, y, z, d, dip
real(8) :: al1, al2, aw1, aw2
real(8) :: disl1, disl2, disl3
! output of DC3D
real(8) :: ux, uy, uz
real(8) :: uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz
integer :: iret
! tmp
real(8) :: dx, dy, t
!real(8) :: alpha ! (lambda+mu)/(lambda+2*mu)
a = (f%lambda+f%mu)/(f%lambda+2.d0*f%mu)
d = f%depth ! 断層上端の深さ。深さは正の値にする。
dip = f%dip
al1 = 0.d0
al2 = f%length ! 断層の長さ。
aw1 = -f%width ! 断層の幅。図を良く見るとAW1はマイナスということで結構紛らわしい。
aw2 = 0.d0
disl1 = f%slip * dcos(f%rake*D2R) ! すべり量を各軸に分解。
disl2 = f%slip * dsin(f%rake*D2R)
disl3 = 0.d0
! 目的の緯度経度をローカル座標系のx, yに変換する。走向の定義が北からなので、ここもややこしい。
call bl2xy(f%lat, f%lon, lat, lon, dx, dy)
t = (-f%str+90.d0)*D2R; ! rotate to local coordinate
x = dx*dcos(t) + dy*dsin(t) ! t = -str+90
y =-dx*dsin(t) + dy*dcos(t)
z =-h
! SUBROUTINE DC3D(ALPHA,X,Y,Z,DEPTH,DIP,
!* AL1,AL2,AW1,AW2,DISL1,DISL2,DISL3,
!* UX,UY,UZ,UXX,UYX,UZX,UXY,UYY,UZY,UXZ,UYZ,UZZ,IRET)
call DC3D(a, x, y, z, d, dip, al1, al2, aw1, aw2, disl1, disl2, disl3, &
ux, uy, uz, uxx, uyx, uzx, uxy, uyy, uzy, uxz, uyz, uzz, iret)
if (iret == 0) then
flt_disp%e = ux*dsin(f%str*D2R) - uy*dcos(f%str*D2R) ! 忘れずに座標変換。走向の向きは北からだからこうなる。
flt_disp%n = ux*dcos(f%str*D2R) + uy*dsin(f%str*D2R)
flt_disp%u = uz
else
! singular location
flt_disp%e = 0.d0
flt_disp%n = 0.d0
flt_disp%u = 0.d0
endif
return
end function flt_disp
! ----------------------------------------------------------
subroutine bl2xy(lat0, lon0, lat, lon, dx, dy)
! ----------------------------------------------------------
! Calculates distance between reference station(lat0, lon0) and
! station2(lat, lon).
!
! Reference: Ohno S.(1987), Geodetic Methods (in Japanese).
! see equations [1.59] and [1.62].
!
! Input:
! lat0, lon0 : position of station1 (reference)
! lat, lon : position of station2
!
! Output:
! dx, dy : distance between stations (sta1 -> sta2).
!
! Created:
! 2013/04/29 Sako
! ----------------------------------------------------------
implicit none
real(8), intent(in) :: lat0, lon0, lat, lon
real(8), intent(out) :: dx, dy
type(ellip) :: ell
real(8) :: dlat, dlon, Bm
real(8) :: a, fi, f, b, e2
! WGS84
ell = ellipsoid('WGS84')
a = ell%a
fi= ell%fi
f = ell%f
b = ell%b
e2= ell%e2
dlat = (lat - lat0)*D2R
dlon = (lon - lon0)*D2R
Bm = dsin(dlat)
dx = a*dcos(lat0*D2R) / dsqrt(1.d0-e2*dsin(lat0*D2R)**2.d0)*dlon
dy = a*(1.d0-e2)*dlat*(1.d0+3.d0/4.d0*e2-3.d0/4.d0*e2*dcos(2.d0*Bm)+1.d0/8.d0*e2*(dlat)**2.d0*dcos(2.d0*Bm))
return
end subroutine bl2xy
! ----------------------------------------------------------
type(ellip) function ellipsoid(ellip_name)
! ----------------------------------------------------------
implicit none
character(len=*), intent(in), optional :: ellip_name
real(8) :: a, fi, f
! default is GRS80
a = 6378137.d0
fi= 298.257222101d0
if (present(ellip_name)) then
if (ellip_name == 'GRS80') then
a = 6378137.d0
fi= 298.257222101d0
elseif (ellip_name == 'WGS84') then
a = 6378137.d0;
fi= 298.257223563d0;
else
! ...
endif
endif
f = 1.d0/fi
ellipsoid%a = a
ellipsoid%fi= fi
ellipsoid%f = f
ellipsoid%b = a*(1.d0-f)
ellipsoid%e2= 1.d0 - (1.d0-f)**2.d0
return
end function ellipsoid
end module okada92_util
2013年4月7日日曜日
fortran90の構造体へのアクセス速度 メモ
前のブログ(http://sak12.blogspot.com/2013/04/matlab.html)で、MATLABの構造体の扱いには注意しなければいけないということをメモしたけれど、これまでfortranで使っていたものは大丈夫なのか不安になった。
で、同じように２乗を計算してみた（1000万回）。
どうやら、fortranでは構造体へのアクセス速度は、特に何も考えなくても大丈夫そう。
構造体をよく使っていたのでドキッとしたのだけれど、一安心。
あと、MATLABが速いことを改めて実感。MATLABでドット演算子を使った場合、fortranのforallと速度はほぼ変わらない。
MATLABは遅いと良く言われたりしてるけど、うまく使えばかなり高速になりそうだし、きちんとMATLABの勉強しようかな。
----- おまけ -----
perlでやってみると、こんなに時間がかかる・・・。まあ数値計算はやっちゃだめということですかね。そりゃそうか。
で、同じように２乗を計算してみた（1000万回）。
どうやら、fortranでは構造体へのアクセス速度は、特に何も考えなくても大丈夫そう。
構造体をよく使っていたのでドキッとしたのだけれど、一安心。
あと、MATLABが速いことを改めて実感。MATLABでドット演算子を使った場合、fortranのforallと速度はほぼ変わらない。
MATLABは遅いと良く言われたりしてるけど、うまく使えばかなり高速になりそうだし、きちんとMATLABの勉強しようかな。
program main
implicit none
type struct1
real(8), allocatable :: x(:)
end type struct1
type struct2
real(8) :: x
end type struct2
integer, parameter :: n = 10000000
type(struct1) :: s1
type(struct2) :: s2(n)
real(8) :: x(n), y(n)
integer :: stime(8), etime(8)
integer :: i
allocate(s1%x(1:n))
forall (i=1:n)
x(i) = 2.d0
s1%x(i) = 2.d0
s2(i)%x = 2.d0
end forall
! x ただの配列同士
forall (i=1:n) y(i) = 0.d0
call date_and_time(values=stime)
forall (i=1:n)
y(i) = x(i) * x(i)
end forall
call date_and_time(values=etime)
write(*,'("elapsed time (x):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0
! s1 構造体が配列を持つ
forall (i=1:n) y(i) = 0.d0
call date_and_time(values=stime)
forall (i=1:n)
y(i) = s1%x(i) * s1%x(i)
end forall
call date_and_time(values=etime)
write(*,'("elapsed time (s1):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0
! s2 構造体配列
forall (i=1:n) y(i) = 0.d0
call date_and_time(values=stime)
forall (i=1:n)
y(i) = s2(i)%x * s2(i)%x
end forall
call date_and_time(values=etime)
write(*,'("elapsed time (s2):", f10.4,"(s)")') dble(etime(7)-stime(7))+ dble(etime(8)-stime(8))/1000.d0
stop
end program main
elapsed time (x): 0.0250(s)
elapsed time (s1): 0.0260(s)
elapsed time (s2): 0.0260(s)
----- おまけ -----
perlでやってみると、こんなに時間がかかる・・・。まあ数値計算はやっちゃだめということですかね。そりゃそうか。
use strict;
use Time::HiRes qw( gettimeofday tv_interval );
my $N = 10_000_000;
my ($i, @x, $y);
foreach $i (1 .. $N) {
push(@x, 2);
}
# estimate square of @x.
my $t0 = [gettimeofday];
for ($i=0;$i<$N;$i++) {
$y = $x[$i]*$x[$i];
}
my ($seconds, $microseconds) = gettimeofday;
my $elapsed = tv_interval($t0, [$seconds, $microseconds]);
print "elapsed time: $elapsed\n";
$ perl test.pl
elapsed time: 1.291174
MATLABのループ メモ
今日は久々にMATLABでコードを書こうと思いたった。
fortran77でよく使われているプログラムがあるのだけど、それをもっと楽に使用するため。
しかし、適当にMATLABでコードを書くと、変なloopを書いてしまいそうで怖い。
で、行列の掛け算のやり方と速度を確認したのでメモしておく。
単純なfor loopと、MATLAB特有のドット演算子を使ったものと、さらに添字にしたらどうなるかな、というもの。
結果は、当然だけどドット演算子を使わないとめちゃくちゃ遅い。あと、添字をつけたらfor loopより遅くなってしまった。使い方がまずいんだろうか。
あと、構造体を絡めて計算する場合は、構造体配列にすると遅くなるので、構造体に配列を入れてしまうのが良いみたい。
fortran77でよく使われているプログラムがあるのだけど、それをもっと楽に使用するため。
しかし、適当にMATLABでコードを書くと、変なloopを書いてしまいそうで怖い。
で、行列の掛け算のやり方と速度を確認したのでメモしておく。
単純なfor loopと、MATLAB特有のドット演算子を使ったものと、さらに添字にしたらどうなるかな、というもの。
結果は、当然だけどドット演算子を使わないとめちゃくちゃ遅い。あと、添字をつけたらfor loopより遅くなってしまった。使い方がまずいんだろうか。
あと、構造体を絡めて計算する場合は、構造体配列にすると遅くなるので、構造体に配列を入れてしまうのが良いみたい。
for loop: y(i) = x(i) * x(i);
Elapsed time is 0.015437 seconds.
vectorized loop: y = x .* x
Elapsed time is 0.002330 seconds.
vectorized loop2: y(1:n) = x(1:n) .* x(1:n)
Elapsed time is 0.021473 seconds.
using struct: y = s.x .* x
Elapsed time is 0.003305 seconds.
using struct: y = s.x .* x
Elapsed time is 0.000005 seconds.
using struct (for loop): y = s2.x .* x
Elapsed time is 0.000239 seconds.
using struct (index): y = s2.x .* x
Elapsed time is 0.001868 seconds.
clear all;
n = 1000000;
s = struct('x', ones(1,n));
x = 2 * ones(1,n);
y = zeros(1,n);
% for loop
disp('for loop: y(i) = x(i) * x(i);');
tic
for i = 1:n
y(i) = x(i) * x(i);
end
toc
% vectorized loop
disp('vectorized loop: y = x .* x');
tic
y = x .* x;
toc
% vectorized loop?
disp('vectorized loop2: y(1:n) = x(1:n) .* x(1:n)');
tic
y(1:n) = x(1:n) .* x(1:n);
toc
% vectorized loop using structure
disp('using struct: y = s.x .* x');
tic
y = s.x .* x;
toc
clear all;
n = 1000;
s = struct('x', ones(1,n)); % s.x(1:n) のようにs.xに配列を指定
s2(1:n) = struct('x', 1); % s(1:n).x のように構造体の配列
y = zeros(1,n);
disp('using struct: y = s.x .* x');
tic
y = s.x .* x;
toc
disp('using struct (for loop): y = s2.x .* x');
tic
s2_x = zeros(1,n);
for i = 1:n
s2_x(i) = s2(i).x;
end
y = s2_x .* x;
toc
disp('using struct (index): y = s2.x .* x');
tic
s2_x = zeros(1,n);
s2_x = s2(1:n).x;
y = s2_x .* x;
toc
登録:
投稿 (Atom)