-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmtrx.f90
140 lines (106 loc) · 4.69 KB
/
mtrx.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
module mtrx
use iso_fortran_env, only: real64
use param, only: CSPEED, N_R, N_THETA, N_PHI, pi, twopi
use epv, only: rp2beta
use dxx, only: cofm
use dmumu, only: cofdu
use file_op, only: open_file_from_environment
use sim3d_utils, only: f0mod
implicit none
contains
! delete this function when the fortran 2008 standards are supported
real(kind=real64) function norm2(arr)
real(kind=real64), intent(in) :: arr(:)
norm2 = sqrt(sum(arr * arr))
end function norm2
function mxptr(gm) result(xptr)
! calculate matrix for xyz' ellipsoid coordinate Kwon to spheric
real(kind=real64), intent(in) :: gm
real(kind=real64) :: xptr(3,3)
real(kind=real64) :: sing, cosg
sing = sin(gm)
cosg = cos(gm)
xptr(1,:) = [0.d0, 0.d0, 1.d0]
xptr(2,:) = [cosg, -sing, 0.d0]
xptr(3,:) = [sing, cosg, 0.d0]
end function
function dmxptr(gm, dgm)
! calculate matrix for xyz' ellipsoid coordinate Kwon to spheric
real(kind=real64), intent(in) :: gm, dgm
real(kind=real64) :: dmxptr(3,3)
real(kind=real64) :: sing, cosg
sing = sin(gm)
cosg = cos(gm)
dmxptr(1,:) = [ 0.d0, 0.d0, 0.d0]
dmxptr(2,:) = [-sing*dgm, -cosg*dgm, 0.d0]
dmxptr(3,:) = [ cosg*dgm, -sing*dgm, 0.d0]
end function
function mbtr(uax1, uax2, uax3) result(b2r)
! calculate martix from magnetic to polar spheric coordinates
real(kind=real64), intent(in) :: uax1(3), uax2(3), uax3(3)
real(kind=real64) :: b2r(3,3)
b2r(:,1) = uax1(:)
b2r(:,2) = uax2(:)
b2r(:,3) = uax3(:)
end function
function mrtx(sintheta, costheta, sinphi, cosphi)
! calculate martix from polar spheric to xyz coordinates
real(kind=real64), intent(in) :: sintheta, costheta, sinphi, cosphi
real(kind=real64) :: mrtx(3,3)
mrtx(1,:) = [sintheta*cosphi, costheta*cosphi, -sinphi]
mrtx(2,:) = [sintheta*sinphi, costheta*sinphi, cosphi]
mrtx(3,:) = [ costheta, -sintheta, 0.d0]
end function
function dmrtx(sintheta, costheta, sinphi, cosphi, dtheta, dphi)
! calculate martix from polar spheric to xyz coordinates
real(kind=real64), intent(in) :: sintheta, costheta, dtheta
real(kind=real64), intent(in) :: sinphi, cosphi, dphi
real(kind=real64) :: dmrtx(3,3)
dmrtx(1,1) = costheta*dtheta*cosphi - sintheta*sinphi*dphi
dmrtx(1,2) = -sintheta*dtheta*cosphi - costheta*sinphi*dphi
dmrtx(1,3) = -cosphi*dphi
dmrtx(2,1) = costheta*dtheta*sinphi + sintheta*cosphi*dphi
dmrtx(2,2) = -sintheta*dtheta*sinphi + costheta*cosphi*dphi
dmrtx(2,3) = -sinphi*dphi
dmrtx(3,1) = -sintheta*dtheta
dmrtx(3,2) = -costheta*dtheta
dmrtx(3,3) = 0.0
end function
! trilinear interpolation
function trilinear(phic, x) result(phi)
! phic the value of phi at the corner of cubic box of side 1
real(kind=real64), intent(in) :: phic(2,2,2)
! location inside the cube (0<=x<=1) or outside x<0 x>1
real(kind=real64), intent(in) :: x(3)
! phi(1)=interpolated phi
real(kind=real64) :: phi
phi = phic(1,1,1) * (1-x(1)) * (1-x(2)) * (1-x(3)) &
+ phic(2,1,1) * x(1) * (1-x(2)) * (1-x(3)) &
+ phic(1,2,1) * (1-x(1)) * x(2) * (1-x(3)) &
+ phic(1,1,2) * (1-x(1)) * (1-x(2)) * x(3) &
+ phic(2,1,2) * x(1) * (1-x(2)) * x(3) &
+ phic(1,2,2) * (1-x(1)) * x(2) * x(3) &
+ phic(2,2,1) * x(1) * x(2) * (1-x(3)) &
+ phic(2,2,2) * x(1) * x(2) * x(3)
end function
function trilineardif(phic, x) result(dphi)
! phic the value of phi at the corner of cubic box of side 1
real(kind=real64), intent(in) :: phic(2,2,2)
! location inside the cube (0<=x<=1) or outside x<0 x>1
real(kind=real64), intent(in) :: x(3)
! dphi(1)=dphi/dx, dphi(2)=dphi/dy, dphi(3)=dphi/dz
real(kind=real64) :: dphi(3)
dphi(1) = (phic(2,1,1)-phic(1,1,1)) * (1-x(2)) * (1-x(3)) &
+ (phic(2,2,1)-phic(1,2,1)) * x(2) * (1-x(3)) &
+ (phic(2,1,2)-phic(1,1,2)) * (1-x(2)) * x(3) &
+ (phic(2,2,2)-phic(1,2,2)) * x(2) * x(3)
dphi(2) = (phic(1,2,1)-phic(1,1,1)) * (1-x(1)) * (1-x(3)) &
+ (phic(2,2,1)-phic(2,1,1)) * x(1) * (1-x(3)) &
+ (phic(1,2,2)-phic(1,1,2)) * (1-x(1)) * x(3) &
+ (phic(2,2,2)-phic(2,1,2)) * x(1) * x(3)
dphi(3) = (phic(1,1,2)-phic(1,1,1)) * (1-x(1)) * (1-x(2)) &
+ (phic(2,1,2)-phic(2,1,1)) * x(1) * (1-x(2)) &
+ (phic(1,2,2)-phic(1,2,1)) * (1-x(1)) * x(2) &
+ (phic(2,2,2)-phic(2,2,1)) * x(1) * x(2)
end function
end module mtrx