-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathdirection_mod.f90
131 lines (64 loc) · 2.8 KB
/
direction_mod.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
module direction_mod
! taken from global_parameters.f90 in ICESHEET
implicit none
double precision, parameter :: pi = 3.141592653589793
contains
double precision function average_direction(direction1, direction2)
implicit none
double precision, intent(in) :: direction1, direction2
if((direction1 > pi/2.0 .and. direction2 < -pi/2.0) ) THEN
average_direction = check_angle( (direction1 + (2.0*pi+direction2))/2.0)
elseif (direction1 < -pi/2.0 .and. direction2 > pi/2.0) THEN
average_direction = check_angle( (direction2 + (2.0*pi+direction1))/2.0)
elseif (direction1 > pi/2.0 .and. (direction2 >-pi/2.0 .and. direction2 < 0.0 )) THEN
if (direction1-direction2 <= pi) THEN
average_direction = (direction1 + direction2) / 2.0
else
average_direction = check_angle( (direction1 + (2.0*pi+direction2))/2.0)
endif
elseif (direction2 > pi/2.0 .and. (direction1 >-pi/2.0 .and. direction1 < 0.0 )) THEN
if (direction2-direction1 <= pi) THEN
average_direction = (direction1 + direction2) / 2.0
else
average_direction = check_angle( (direction2 + (2.0*pi+direction1))/2.0)
endif
else
average_direction = (direction1 + direction2) / 2.0
endif
end function average_direction
double precision function check_angle(angle)
implicit none
double precision, intent(in) :: angle
! double precision :: check_angle
if (angle > pi) THEN
check_angle = angle - 2.0 * pi
elseif (angle < -pi) THEN
check_angle = angle + 2.0 * pi
else
check_angle = angle
end if
end function check_angle
double precision function general_direction(direction1, direction2, rotate_factor)
implicit none
double precision, intent(in) :: direction1, direction2, rotate_factor
if((direction1 >= pi/2.0 .and. direction2 <= -pi/2.0) ) THEN
general_direction = check_angle( ((2*pi+direction2)-direction1)*rotate_factor + direction1)
elseif (direction1 <= -pi/2.0 .and. direction2 >= pi/2.0) THEN
general_direction = check_angle( direction1 - ((2.0*pi+direction1)-direction2)*rotate_factor )
elseif (direction1 >= pi/2.0 .and. (direction2 >=-pi/2.0 .and. direction2 < 0.0 )) THEN
if (direction1-direction2 <= pi) THEN
general_direction = direction1 - (direction1 - direction2) *rotate_factor
else
general_direction = check_angle( direction1+(2.0*pi -(direction1-direction2))*rotate_factor)
endif
elseif (direction2 >= pi/2.0 .and. (direction1 >=-pi/2.0 .and. direction1 <= 0.0 )) THEN
if (direction2-direction1 <= pi) THEN
general_direction = direction1 + (direction2 - direction1)*rotate_factor
else
general_direction = check_angle(direction1 - (2.0*pi -(direction2 - direction1))*rotate_factor)
endif
else
general_direction = check_angle(direction1 - (direction1-direction2)*rotate_factor)
endif
end function general_direction
end module direction_mod