| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | module raffle__misc_linalg | ||
| 2 | !! Module contains various linear algebra functions and subroutines. | ||
| 3 | use raffle__constants, only: real32, pi | ||
| 4 | implicit none | ||
| 5 | |||
| 6 | |||
| 7 | private | ||
| 8 | |||
| 9 | public :: modu, cross | ||
| 10 | public :: get_distance, get_angle, get_dihedral_angle | ||
| 11 | public :: get_improper_dihedral_angle | ||
| 12 | public :: inverse_3x3 | ||
| 13 | |||
| 14 | |||
| 15 | interface get_angle | ||
| 16 | procedure get_angle_from_points, get_angle_from_vectors | ||
| 17 | end interface get_angle | ||
| 18 | |||
| 19 | interface get_dihedral_angle | ||
| 20 | procedure get_dihedral_angle_from_points, get_dihedral_angle_from_vectors | ||
| 21 | end interface get_dihedral_angle | ||
| 22 | |||
| 23 | interface get_improper_dihedral_angle | ||
| 24 | procedure get_improper_dihedral_angle_from_points, & | ||
| 25 | get_improper_dihedral_angle_from_vectors | ||
| 26 | end interface get_improper_dihedral_angle | ||
| 27 | |||
| 28 | contains | ||
| 29 | |||
| 30 | !############################################################################### | ||
| 31 |
1/2✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
|
14 | pure function modu(vector) |
| 32 | !! Return the magnitude of a vector of any size. | ||
| 33 | implicit none | ||
| 34 | |||
| 35 | ! Arguments | ||
| 36 | real(real32),dimension(:), intent(in) :: vector | ||
| 37 | !! Input vector. | ||
| 38 | real(real32) :: modu | ||
| 39 | !! Output magnitude. | ||
| 40 | |||
| 41 |
2/2✓ Branch 0 taken 42 times.
✓ Branch 1 taken 14 times.
|
56 | modu = sqrt( sum( vector ** 2 ) ) |
| 42 | 14 | end function modu | |
| 43 | !############################################################################### | ||
| 44 | |||
| 45 | |||
| 46 | !############################################################################### | ||
| 47 |
1/2✓ Branch 0 taken 434215091 times.
✗ Branch 1 not taken.
|
434215091 | pure function cross(a,b) |
| 48 | !! Return the cross product of two vectors. | ||
| 49 | implicit none | ||
| 50 | |||
| 51 | ! Arguments | ||
| 52 | real(real32), dimension(3), intent(in) :: a, b | ||
| 53 | !! Input vectors. | ||
| 54 | real(real32), dimension(3) :: cross | ||
| 55 | !! Output cross product. | ||
| 56 | |||
| 57 | 434215091 | cross(1) = a(2) * b(3) - a(3) * b(2) | |
| 58 | 434215091 | cross(2) = a(3) * b(1) - a(1) * b(3) | |
| 59 | 434215091 | cross(3) = a(1) * b(2) - a(2) * b(1) | |
| 60 | |||
| 61 | end function cross | ||
| 62 | !############################################################################### | ||
| 63 | |||
| 64 | |||
| 65 | !############################################################################### | ||
| 66 | 1 | pure function get_distance(point1, point2) result(distance) | |
| 67 | !! Return the distance between two points. | ||
| 68 | implicit none | ||
| 69 | |||
| 70 | ! Arguments | ||
| 71 | real(real32), dimension(3), intent(in) :: point1, point2 | ||
| 72 | !! Input points. | ||
| 73 | real(real32) :: distance | ||
| 74 | !! Output distance. | ||
| 75 | |||
| 76 |
4/6✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
4 | distance = norm2( point1 - point2 ) |
| 77 | |||
| 78 | 1 | return | |
| 79 | end function get_distance | ||
| 80 | !############################################################################### | ||
| 81 | |||
| 82 | |||
| 83 | !############################################################################### | ||
| 84 | 216879957 | pure function get_angle_from_vectors(vector1, vector2) result(angle) | |
| 85 | !! Return the angle between two vectors. | ||
| 86 | implicit none | ||
| 87 | |||
| 88 | ! Arguments | ||
| 89 | real(real32), dimension(3), intent(in) :: vector1, vector2 | ||
| 90 | !! Input vectors. | ||
| 91 | real(real32) :: angle | ||
| 92 | !! Output angle. | ||
| 93 | |||
| 94 | angle = dot_product(vector1,vector2) / & | ||
| 95 |
14/14✓ Branch 0 taken 650639871 times.
✓ Branch 1 taken 216879957 times.
✓ Branch 2 taken 650639871 times.
✓ Branch 3 taken 216879957 times.
✓ Branch 4 taken 607240982 times.
✓ Branch 5 taken 43398889 times.
✓ Branch 6 taken 267338152 times.
✓ Branch 7 taken 339902830 times.
✓ Branch 8 taken 650639871 times.
✓ Branch 9 taken 216879957 times.
✓ Branch 10 taken 620032204 times.
✓ Branch 11 taken 30607667 times.
✓ Branch 12 taken 328552269 times.
✓ Branch 13 taken 291479935 times.
|
2168799570 | ( norm2(vector1) * norm2(vector2) ) |
| 96 |
2/2✓ Branch 0 taken 2120909 times.
✓ Branch 1 taken 214759048 times.
|
216879957 | if(angle .ge. 1._real32)then |
| 97 | 2120909 | angle = 0._real32 | |
| 98 |
2/2✓ Branch 0 taken 2111717 times.
✓ Branch 1 taken 212647331 times.
|
214759048 | elseif(angle .le. -1._real32)then |
| 99 | 2111717 | angle = pi | |
| 100 | else | ||
| 101 | 212647331 | angle = acos(angle) | |
| 102 | end if | ||
| 103 | 216879957 | end function get_angle_from_vectors | |
| 104 | !############################################################################### | ||
| 105 | |||
| 106 | |||
| 107 | !############################################################################### | ||
| 108 | 5694883 | pure function get_angle_from_points(point1, point2, point3) result(angle) | |
| 109 | !! Return the angle formed by three points. | ||
| 110 | !! | ||
| 111 | !! The angle is formed by the path point1 -> point2 -> point3. | ||
| 112 | implicit none | ||
| 113 | |||
| 114 | ! Arguments | ||
| 115 | real(real32), dimension(3), intent(in) :: point1, point2, point3 | ||
| 116 | !! Input points. | ||
| 117 | real(real32) :: angle | ||
| 118 | !! Output angle. | ||
| 119 | |||
| 120 | angle = dot_product( point1 - point2, point3 - point2 ) / & | ||
| 121 |
14/14✓ Branch 0 taken 17084649 times.
✓ Branch 1 taken 5694883 times.
✓ Branch 2 taken 17084649 times.
✓ Branch 3 taken 5694883 times.
✓ Branch 4 taken 16514465 times.
✓ Branch 5 taken 570184 times.
✓ Branch 6 taken 7352792 times.
✓ Branch 7 taken 9161673 times.
✓ Branch 8 taken 17084649 times.
✓ Branch 9 taken 5694883 times.
✓ Branch 10 taken 16502117 times.
✓ Branch 11 taken 582532 times.
✓ Branch 12 taken 6890367 times.
✓ Branch 13 taken 9611750 times.
|
56948830 | ( norm2( point1 - point2 ) * norm2( point3 - point2 ) ) |
| 122 |
2/2✓ Branch 0 taken 20 times.
✓ Branch 1 taken 5694863 times.
|
5694883 | if(angle .ge. 1._real32)then |
| 123 | 20 | angle = 0._real32 | |
| 124 |
2/2✓ Branch 0 taken 2698 times.
✓ Branch 1 taken 5692165 times.
|
5694863 | elseif(angle .le. -1._real32)then |
| 125 | 2698 | angle = pi | |
| 126 | else | ||
| 127 | 5692165 | angle = acos(angle) | |
| 128 | end if | ||
| 129 | 5694883 | end function get_angle_from_points | |
| 130 | !############################################################################### | ||
| 131 | |||
| 132 | |||
| 133 | !############################################################################### | ||
| 134 | 1 | pure function get_dihedral_angle_from_vectors( & | |
| 135 | vector1, vector2, vector3) result(angle) | ||
| 136 | !! Return the dihedral angle between two planes. | ||
| 137 | !! | ||
| 138 | !! The dihedral angle is the angle between the plane defined by the cross | ||
| 139 | !! product of two vectors and a third vector. | ||
| 140 | !! i.e. ( vector1 x vector2 ) . vector3 | ||
| 141 | implicit none | ||
| 142 | |||
| 143 | ! Arguments | ||
| 144 | real(real32), dimension(3), intent(in) :: vector1, vector2, vector3 | ||
| 145 | !! Input vectors. | ||
| 146 | real(real32) :: angle | ||
| 147 | !! Output angle. | ||
| 148 | |||
| 149 | 1 | angle = get_angle(cross(vector1, vector2), vector3) | |
| 150 | |||
| 151 | 1 | end function get_dihedral_angle_from_vectors | |
| 152 | !############################################################################### | ||
| 153 | |||
| 154 | |||
| 155 | !############################################################################### | ||
| 156 | 1 | pure function get_dihedral_angle_from_points( & | |
| 157 | point1, point2, point3, point4 & | ||
| 158 | ) result(angle) | ||
| 159 | !! Return the dihedral angle between two planes. | ||
| 160 | !! | ||
| 161 | !! The dihedral angle is the angle between the plane defined by four points. | ||
| 162 | !! i.e. ( point2 - point1 ) x ( point3 - point2 ) . ( point4 - point2 ) | ||
| 163 | !! alt. angle between plane point1point2point3 and vector point2point4 | ||
| 164 | implicit none | ||
| 165 | real(real32), dimension(3), intent(in) :: point1, point2, point3, point4 | ||
| 166 | real(real32) :: angle | ||
| 167 | |||
| 168 |
6/6✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
✓ Branch 3 taken 1 times.
✓ Branch 5 taken 3 times.
✓ Branch 6 taken 1 times.
|
10 | angle = get_angle(cross(point2 - point1, point3 - point2), point4 - point2) |
| 169 | |||
| 170 | 1 | end function get_dihedral_angle_from_points | |
| 171 | !############################################################################### | ||
| 172 | |||
| 173 | |||
| 174 | !############################################################################### | ||
| 175 | 478782 | pure function get_improper_dihedral_angle_from_vectors( & | |
| 176 | vector1, vector2, vector3 & | ||
| 177 | ) result(angle) | ||
| 178 | !! Return the improper dihedral angle between two planes. | ||
| 179 | !! | ||
| 180 | !! The improper dihedral angle is the angle between two planes made by | ||
| 181 | !! three vectors. | ||
| 182 | !! i.e. ( vector1 x vector2 ) . ( vector2 x vector3 ) | ||
| 183 | !! alt. angle between plane vector1vector2 and vector2vector3 | ||
| 184 | implicit none | ||
| 185 | real(real32), dimension(3), intent(in) :: vector1, vector2, vector3 | ||
| 186 | real(real32) :: angle | ||
| 187 | |||
| 188 | angle = get_angle( & | ||
| 189 | cross(vector1, vector2), & | ||
| 190 | cross(vector2, vector3) & | ||
| 191 | 478782 | ) | |
| 192 | !! map angle back into the range [0, pi] | ||
| 193 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 478782 times.
|
478782 | if(angle .gt. pi) angle = 2._real32 * pi - angle |
| 194 | |||
| 195 | |||
| 196 | 478782 | end function get_improper_dihedral_angle_from_vectors | |
| 197 | !############################################################################### | ||
| 198 | |||
| 199 | |||
| 200 | !############################################################################### | ||
| 201 | 216393559 | pure function get_improper_dihedral_angle_from_points( & | |
| 202 | point1, point2, point3, point4 & | ||
| 203 | ) result(angle) | ||
| 204 | !! Return the improper dihedral angle between two planes. | ||
| 205 | !! | ||
| 206 | !! The dihedral angle is the angle between the plane defined by four points. | ||
| 207 | !! i.e. ( point2 - point1 ) x ( point3 - point1 ) . | ||
| 208 | !! ( point4 - point2 ) x ( point3 - point1 ) | ||
| 209 | !! alt. angle between plane point1point2point3 and point1point3point4 | ||
| 210 | implicit none | ||
| 211 | real(real32), dimension(3), intent(in) :: point1, point2, point3, point4 | ||
| 212 | real(real32) :: angle | ||
| 213 | |||
| 214 | angle = get_angle( & | ||
| 215 | cross(point2 - point1, point3 - point1), & | ||
| 216 | cross(point3 - point1, point4 - point1) & | ||
| 217 |
8/8✓ Branch 0 taken 649180677 times.
✓ Branch 1 taken 216393559 times.
✓ Branch 2 taken 649180677 times.
✓ Branch 3 taken 216393559 times.
✓ Branch 5 taken 649180677 times.
✓ Branch 6 taken 216393559 times.
✓ Branch 7 taken 649180677 times.
✓ Branch 8 taken 216393559 times.
|
2813116267 | ) |
| 218 | |||
| 219 | 216393559 | end function get_improper_dihedral_angle_from_points | |
| 220 | !############################################################################### | ||
| 221 | |||
| 222 | |||
| 223 | !############################################################################### | ||
| 224 |
1/2✓ Branch 0 taken 78411 times.
✗ Branch 1 not taken.
|
78411 | pure function inverse_3x3(mat) result(output) |
| 225 | implicit none | ||
| 226 | real(real32) :: det | ||
| 227 | real(real32), dimension(3,3) :: output | ||
| 228 | real(real32), dimension(3,3), intent(in) :: mat | ||
| 229 | |||
| 230 | det = & | ||
| 231 | mat(1,1) * mat(2,2) * mat(3,3) - mat(1,1) * mat(2,3) * mat(3,2) - & | ||
| 232 | mat(1,2) * mat(2,1) * mat(3,3) + mat(1,2) * mat(2,3) * mat(3,1) + & | ||
| 233 | 78411 | mat(1,3) * mat(2,1) * mat(3,2) - mat(1,3) * mat(2,2) * mat(3,1) | |
| 234 | |||
| 235 | 78411 | output(1,1) = +1._real32 / det * (mat(2,2) * mat(3,3) - mat(2,3) * mat(3,2)) | |
| 236 | 78411 | output(2,1) = -1._real32 / det * (mat(2,1) * mat(3,3) - mat(2,3) * mat(3,1)) | |
| 237 | 78411 | output(3,1) = +1._real32 / det * (mat(2,1) * mat(3,2) - mat(2,2) * mat(3,1)) | |
| 238 | 78411 | output(1,2) = -1._real32 / det * (mat(1,2) * mat(3,3) - mat(1,3) * mat(3,2)) | |
| 239 | 78411 | output(2,2) = +1._real32 / det * (mat(1,1) * mat(3,3) - mat(1,3) * mat(3,1)) | |
| 240 | 78411 | output(3,2) = -1._real32 / det * (mat(1,1) * mat(3,2) - mat(1,2) * mat(3,1)) | |
| 241 | 78411 | output(1,3) = +1._real32 / det * (mat(1,2) * mat(2,3) - mat(1,3) * mat(2,2)) | |
| 242 | 78411 | output(2,3) = -1._real32 / det * (mat(1,1) * mat(2,3) - mat(1,3) * mat(2,1)) | |
| 243 | 78411 | output(3,3) = +1._real32 / det * (mat(1,1) * mat(2,2) - mat(1,2) * mat(2,1)) | |
| 244 | |||
| 245 | 78411 | end function inverse_3x3 | |
| 246 | !############################################################################### | ||
| 247 | |||
| 248 | end module raffle__misc_linalg | ||
| 249 |