DAMASK with grid solvers  Revision: v2.0.3-2204-gdb1f2151
The Düsseldorf Advanced Material Simulation Kit with Grid Solvers
rotations.f90
Go to the documentation of this file.
1 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/rotations.f90"
2 # 1 "<built-in>"
3 # 1 "<command-line>"
4 # 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/rotations.f90"
5 ! ###################################################################
6 ! Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
7 ! Modified 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
8 ! All rights reserved.
9 !
10 ! Redistribution and use in source and binary forms, with or without modification, are
11 ! permitted provided that the following conditions are met:
12 !
13 ! - Redistributions of source code must retain the above copyright notice, this list
14 ! of conditions and the following disclaimer.
15 ! - Redistributions in binary form must reproduce the above copyright notice, this
16 ! list of conditions and the following disclaimer in the documentation and/or
17 ! other materials provided with the distribution.
18 ! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
19 ! of its contributors may be used to endorse or promote products derived from
20 ! this software without specific prior written permission.
21 !
22 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
23 ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24 ! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25 ! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
26 ! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
27 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
28 ! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
29 ! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
30 ! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
31 ! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
32 ! ###################################################################
33 
34 !---------------------------------------------------------------------------------------------------
40 !
41 ! All methods and naming conventions based on Rowenhorst_etal2015
42 ! Convention 1: coordinate frames are right-handed
43 ! Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
44 ! when viewing from the end point of the rotation axis towards the origin
45 ! Convention 3: rotations will be interpreted in the passive sense
46 ! Convention 4: Euler angle triplets are implemented using the Bunge convention,
47 ! with the angular ranges as [0, 2π],[0, π],[0, 2π]
48 ! Convention 5: the rotation angle ω is limited to the interval [0, π]
49 ! Convention 6: the real part of a quaternion is positive, Re(q) > 0
50 ! Convention 7: P = -1
51 !---------------------------------------------------------------------------------------------------
52 
53 module rotations
54  use prec
55  use io
56  use math
57  use lambert
58  use quaternions
59 
60  implicit none
61  private
62 
63  type, public :: rotation
64  type(quaternion), private :: q
65  contains
66  procedure, public :: asquaternion
67  procedure, public :: aseulers
68  procedure, public :: asaxisangle
69  procedure, public :: asrodrigues
70  procedure, public :: asmatrix
71  !------------------------------------------
72  procedure, public :: fromquaternion
73  procedure, public :: fromeulers
74  procedure, public :: fromaxisangle
75  procedure, public :: frommatrix
76  !------------------------------------------
77  procedure, private :: rotrot__
78  generic, public :: operator(*) => rotrot__
79  generic, public :: rotate => rotvector,rottensor2,rottensor4
80  procedure, public :: rotvector
81  procedure, public :: rottensor2
82  procedure, public :: rottensor4
83  procedure, public :: rottensor4sym
84  procedure, public :: misorientation
85  procedure, public :: standardize
86  end type rotation
87 
88  public :: &
89  rotations_init, &
90  eu2om
91 
92 contains
93 
94 !--------------------------------------------------------------------------------------------------
96 !--------------------------------------------------------------------------------------------------
97 subroutine rotations_init
98 
99  call quaternions_init
100  write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6)
101  call unittest
102 
103 end subroutine rotations_init
104 
105 
106 !---------------------------------------------------------------------------------------------------
107 ! Return rotation in different representations
108 !---------------------------------------------------------------------------------------------------
109 pure function asquaternion(self)
110 
111  class(rotation), intent(in) :: self
112  real(preal), dimension(4) :: asquaternion
113 
114  asquaternion = self%q%asArray()
115 
116 end function asquaternion
117 !---------------------------------------------------------------------------------------------------
118 pure function aseulers(self)
119 
120  class(rotation), intent(in) :: self
121  real(preal), dimension(3) :: aseulers
122 
123  aseulers = qu2eu(self%q%asArray())
124 
125 end function aseulers
126 !---------------------------------------------------------------------------------------------------
127 pure function asaxisangle(self)
128 
129  class(rotation), intent(in) :: self
130  real(preal), dimension(4) :: asaxisangle
131 
132  asaxisangle = qu2ax(self%q%asArray())
133 
134 end function asaxisangle
135 !---------------------------------------------------------------------------------------------------
136 pure function asmatrix(self)
137 
138  class(rotation), intent(in) :: self
139  real(preal), dimension(3,3) :: asmatrix
140 
141  asmatrix = qu2om(self%q%asArray())
142 
143 end function asmatrix
144 !---------------------------------------------------------------------------------------------------
145 pure function asrodrigues(self)
146 
147  class(rotation), intent(in) :: self
148  real(preal), dimension(4) :: asrodrigues
149 
150  asrodrigues = qu2ro(self%q%asArray())
151 
152 end function asrodrigues
153 !---------------------------------------------------------------------------------------------------
154 pure function ashomochoric(self)
155 
156  class(rotation), intent(in) :: self
157  real(preal), dimension(3) :: ashomochoric
158 
159  ashomochoric = qu2ho(self%q%asArray())
160 
161 end function ashomochoric
162 
163 !---------------------------------------------------------------------------------------------------
164 ! Initialize rotation from different representations
165 !---------------------------------------------------------------------------------------------------
166 subroutine fromquaternion(self,qu)
167 
168  class(rotation), intent(out) :: self
169  real(pReal), dimension(4), intent(in) :: qu
170 
171  if (dneq(norm2(qu),1.0_preal)) &
172  call io_error(402,ext_msg='fromQuaternion')
173 
174  self%q = qu
175 
176 end subroutine fromquaternion
177 !---------------------------------------------------------------------------------------------------
178 subroutine fromeulers(self,eu,degrees)
179 
180  class(rotation), intent(out) :: self
181  real(pReal), dimension(3), intent(in) :: eu
182  logical, intent(in), optional :: degrees
183 
184  real(pReal), dimension(3) :: Eulers
185 
186  if (.not. present(degrees)) then
187  eulers = eu
188  else
189  eulers = merge(eu*inrad,eu,degrees)
190  endif
191 
192  if (any(eulers<0.0_preal) .or. any(eulers>2.0_preal*pi) .or. eulers(2) > pi) &
193  call io_error(402,ext_msg='fromEulers')
194 
195  self%q = eu2qu(eulers)
196 
197 end subroutine fromeulers
198 !---------------------------------------------------------------------------------------------------
199 subroutine fromaxisangle(self,ax,degrees,P)
200 
201  class(rotation), intent(out) :: self
202  real(pReal), dimension(4), intent(in) :: ax
203  logical, intent(in), optional :: degrees
204  integer, intent(in), optional :: P
205 
206  real(pReal) :: angle
207  real(pReal),dimension(3) :: axis
208 
209  if (.not. present(degrees)) then
210  angle = ax(4)
211  else
212  angle = merge(ax(4)*inrad,ax(4),degrees)
213  endif
214 
215  if (.not. present(p)) then
216  axis = ax(1:3)
217  else
218  axis = ax(1:3) * merge(-1.0_preal,1.0_preal,p == 1)
219  if(abs(p) /= 1) call io_error(402,ext_msg='fromAxisAngle (P)')
220  endif
221 
222  if(dneq(norm2(axis),1.0_preal) .or. angle < 0.0_preal .or. angle > pi) &
223  call io_error(402,ext_msg='fromAxisAngle')
224 
225  self%q = ax2qu([axis,angle])
226 
227 end subroutine fromaxisangle
228 !---------------------------------------------------------------------------------------------------
229 subroutine frommatrix(self,om)
230 
231  class(rotation), intent(out) :: self
232  real(pReal), dimension(3,3), intent(in) :: om
233 
234  if (dneq(math_det33(om),1.0_preal,tol=1.0e-5_preal)) &
235  call io_error(402,ext_msg='fromMatrix')
236 
237  self%q = om2qu(om)
238 
239 end subroutine frommatrix
240 !---------------------------------------------------------------------------------------------------
241 
242 
243 !---------------------------------------------------------------------------------------------------
245 !---------------------------------------------------------------------------------------------------
246 pure elemental function rotrot__(self,R) result(rRot)
247 
248  type(rotation) :: rrot
249  class(rotation), intent(in) :: self,r
250 
251  rrot = rotation(self%q*r%q)
252  call rrot%standardize()
253 
254 end function rotrot__
255 
256 
257 !---------------------------------------------------------------------------------------------------
259 !---------------------------------------------------------------------------------------------------
260 pure elemental subroutine standardize(self)
261 
262  class(rotation), intent(inout) :: self
263 
264  if (real(self%q) < 0.0_preal) self%q = self%q%homomorphed()
265 
266 end subroutine standardize
267 
268 
269 !---------------------------------------------------------------------------------------------------
272 !---------------------------------------------------------------------------------------------------
273 pure function rotvector(self,v,active) result(vRot)
274 
275  real(preal), dimension(3) :: vrot
276  class(rotation), intent(in) :: self
277  real(preal), intent(in), dimension(3) :: v
278  logical, intent(in), optional :: active
279 
280  real(preal), dimension(3) :: v_normed
281  type(quaternion) :: q
282  logical :: passive
283 
284  if (present(active)) then
285  passive = .not. active
286  else
287  passive = .true.
288  endif
289 
290  if (deq0(norm2(v))) then
291  vrot = v
292  else
293  v_normed = v/norm2(v)
294  if (passive) then
295  q = self%q * (quaternion([0.0_preal, v_normed(1), v_normed(2), v_normed(3)]) * conjg(self%q) )
296  else
297  q = conjg(self%q) * (quaternion([0.0_preal, v_normed(1), v_normed(2), v_normed(3)]) * self%q )
298  endif
299  vrot = q%aimag()*norm2(v)
300  endif
301 
302 end function rotvector
303 
304 
305 !---------------------------------------------------------------------------------------------------
309 !---------------------------------------------------------------------------------------------------
310 pure function rottensor2(self,T,active) result(tRot)
311 
312  real(preal), dimension(3,3) :: trot
313  class(rotation), intent(in) :: self
314  real(preal), intent(in), dimension(3,3) :: t
315  logical, intent(in), optional :: active
316 
317  logical :: passive
318 
319  if (present(active)) then
320  passive = .not. active
321  else
322  passive = .true.
323  endif
324 
325  if (passive) then
326  trot = matmul(matmul(self%asMatrix(),t),transpose(self%asMatrix()))
327  else
328  trot = matmul(matmul(transpose(self%asMatrix()),t),self%asMatrix())
329  endif
330 
331 end function rottensor2
332 
333 
334 !---------------------------------------------------------------------------------------------------
339 !---------------------------------------------------------------------------------------------------
340 pure function rottensor4(self,T,active) result(tRot)
341 
342  real(preal), dimension(3,3,3,3) :: trot
343  class(rotation), intent(in) :: self
344  real(preal), intent(in), dimension(3,3,3,3) :: t
345  logical, intent(in), optional :: active
346 
347  real(preal), dimension(3,3) :: r
348  integer :: i,j,k,l,m,n,o,p
349 
350  if (present(active)) then
351  r = merge(transpose(self%asMatrix()),self%asMatrix(),active)
352  else
353  r = self%asMatrix()
354  endif
355 
356  trot = 0.0_preal
357  do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
358  do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
359  trot(i,j,k,l) = trot(i,j,k,l) &
360  + r(i,m) * r(j,n) * r(k,o) * r(l,p) * t(m,n,o,p)
361  enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
362 
363 end function rottensor4
364 
365 
366 !---------------------------------------------------------------------------------------------------
370 !---------------------------------------------------------------------------------------------------
371 pure function rottensor4sym(self,T,active) result(tRot)
372 
373  real(preal), dimension(6,6) :: trot
374  class(rotation), intent(in) :: self
375  real(preal), intent(in), dimension(6,6) :: t
376  logical, intent(in), optional :: active
377 
378  if (present(active)) then
379  trot = math_sym3333to66(rottensor4(self,math_66tosym3333(t),active))
380  else
381  trot = math_sym3333to66(rottensor4(self,math_66tosym3333(t)))
382  endif
383 
384 end function rottensor4sym
385 
386 
387 !---------------------------------------------------------------------------------------------------
389 !---------------------------------------------------------------------------------------------------
390 pure elemental function misorientation(self,other)
391 
392  type(rotation) :: misorientation
393  class(rotation), intent(in) :: self, other
394 
395  misorientation%q = other%q * conjg(self%q)
396 
397 end function misorientation
398 
399 
400 !---------------------------------------------------------------------------------------------------
403 !---------------------------------------------------------------------------------------------------
404 pure function qu2om(qu) result(om)
405 
406  real(preal), intent(in), dimension(4) :: qu
407  real(preal), dimension(3,3) :: om
408 
409  real(preal) :: qq
410 
411  qq = qu(1)**2-sum(qu(2:4)**2)
412 
413 
414  om(1,1) = qq+2.0_preal*qu(2)**2
415  om(2,2) = qq+2.0_preal*qu(3)**2
416  om(3,3) = qq+2.0_preal*qu(4)**2
417 
418  om(1,2) = 2.0_preal*(qu(2)*qu(3)-qu(1)*qu(4))
419  om(2,3) = 2.0_preal*(qu(3)*qu(4)-qu(1)*qu(2))
420  om(3,1) = 2.0_preal*(qu(4)*qu(2)-qu(1)*qu(3))
421  om(2,1) = 2.0_preal*(qu(3)*qu(2)+qu(1)*qu(4))
422  om(3,2) = 2.0_preal*(qu(4)*qu(3)+qu(1)*qu(2))
423  om(1,3) = 2.0_preal*(qu(2)*qu(4)+qu(1)*qu(3))
424 
425  if (p < 0.0_preal) om = transpose(om)
426 
427 end function qu2om
428 
429 
430 !---------------------------------------------------------------------------------------------------
433 !---------------------------------------------------------------------------------------------------
434 pure function qu2eu(qu) result(eu)
435 
436  real(preal), intent(in), dimension(4) :: qu
437  real(preal), dimension(3) :: eu
438 
439  real(preal) :: q12, q03, chi, chiinv
440 
441  q03 = qu(1)**2+qu(4)**2
442  q12 = qu(2)**2+qu(3)**2
443  chi = sqrt(q03*q12)
444 
445  degenerated: if (deq0(chi)) then
446  eu = merge([atan2(-p*2.0_preal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_preal, 0.0_preal], &
447  [atan2( 2.0_preal*qu(2)*qu(3),qu(2)**2-qu(3)**2), pi, 0.0_preal], &
448  deq0(q12))
449  else degenerated
450  chiinv = 1.0_preal/chi
451  eu = [atan2((-p*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-p*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
452  atan2( 2.0_preal*chi, q03-q12 ), &
453  atan2(( p*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-p*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
454  endif degenerated
455  where(eu<0.0_preal) eu = mod(eu+2.0_preal*pi,[2.0_preal*pi,pi,2.0_preal*pi])
456 
457 end function qu2eu
458 
459 
460 !---------------------------------------------------------------------------------------------------
463 !---------------------------------------------------------------------------------------------------
464 pure function qu2ax(qu) result(ax)
465 
466  real(preal), intent(in), dimension(4) :: qu
467  real(preal), dimension(4) :: ax
468 
469  real(preal) :: omega, s
470 
471  if (deq0(sum(qu(2:4)**2))) then
472  ax = [ 0.0_preal, 0.0_preal, 1.0_preal, 0.0_preal ] ! axis = [001]
473  elseif (dneq0(qu(1))) then
474  s = sign(1.0_preal,qu(1))/norm2(qu(2:4))
475  omega = 2.0_preal * acos(math_clip(qu(1),-1.0_preal,1.0_preal))
476  ax = [ qu(2)*s, qu(3)*s, qu(4)*s, omega ]
477  else
478  ax = [ qu(2), qu(3), qu(4), pi ]
479  end if
480 
481 end function qu2ax
482 
483 
484 !---------------------------------------------------------------------------------------------------
487 !---------------------------------------------------------------------------------------------------
488 pure function qu2ro(qu) result(ro)
489 
490  real(preal), intent(in), dimension(4) :: qu
491  real(preal), dimension(4) :: ro
492 
493  real(preal) :: s
494  real(preal), parameter :: thr = 1.0e-8_preal
495 
496  if (abs(qu(1)) < thr) then
497  ro = [qu(2), qu(3), qu(4), ieee_value(1.0_preal,ieee_positive_inf)]
498  else
499  s = norm2(qu(2:4))
500  if (s < thr) then
501  ro = [0.0_preal, 0.0_preal, p, 0.0_preal]
502  else
503  ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_preal,1.0_preal)))]
504  endif
505 
506  end if
507 
508 end function qu2ro
509 
510 
511 !---------------------------------------------------------------------------------------------------
514 !---------------------------------------------------------------------------------------------------
515 pure function qu2ho(qu) result(ho)
516 
517  real(preal), intent(in), dimension(4) :: qu
518  real(preal), dimension(3) :: ho
519 
520  real(preal) :: omega, f
521 
522  omega = 2.0 * acos(math_clip(qu(1),-1.0_preal,1.0_preal))
523 
524  if (deq0(omega)) then
525  ho = [ 0.0_preal, 0.0_preal, 0.0_preal ]
526  else
527  ho = qu(2:4)
528  f = 0.75_preal * ( omega - sin(omega) )
529  ho = ho/norm2(ho)* f**(1.0_preal/3.0_preal)
530  end if
531 
532 end function qu2ho
533 
534 
535 !---------------------------------------------------------------------------------------------------
538 !---------------------------------------------------------------------------------------------------
539 pure function qu2cu(qu) result(cu)
540 
541  real(preal), intent(in), dimension(4) :: qu
542  real(preal), dimension(3) :: cu
543 
544  cu = ho2cu(qu2ho(qu))
545 
546 end function qu2cu
547 
548 
549 !---------------------------------------------------------------------------------------------------
553 !---------------------------------------------------------------------------------------------------
554 pure function om2qu(om) result(qu)
555 
556  real(preal), intent(in), dimension(3,3) :: om
557  real(preal), dimension(4) :: qu
558 
559  qu = eu2qu(om2eu(om))
560 
561 end function om2qu
562 
563 
564 !---------------------------------------------------------------------------------------------------
567 !---------------------------------------------------------------------------------------------------
568 pure function om2eu(om) result(eu)
569 
570  real(preal), intent(in), dimension(3,3) :: om
571  real(preal), dimension(3) :: eu
572  real(preal) :: zeta
573 
574  if (abs(om(3,3)) < 1.0_preal) then
575  zeta = 1.0_preal/sqrt(1.0_preal-om(3,3)**2.0_preal)
576  eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), &
577  acos(om(3,3)), &
578  atan2(om(1,3)*zeta, om(2,3)*zeta)]
579  else
580  eu = [atan2(om(1,2),om(1,1)), 0.5_preal*pi*(1.0_preal-om(3,3)),0.0_preal ]
581  end if
582 
583  where(eu<0.0_preal) eu = mod(eu+2.0_preal*pi,[2.0_preal*pi,pi,2.0_preal*pi])
584 
585 end function om2eu
586 
587 
588 !---------------------------------------------------------------------------------------------------
591 !---------------------------------------------------------------------------------------------------
592 function om2ax(om) result(ax)
593 
594  real(preal), intent(in), dimension(3,3) :: om
595  real(preal), dimension(4) :: ax
596 
597  real(preal) :: t
598  real(preal), dimension(3) :: wr, wi
599  real(preal), dimension((64+2)*3) :: work
600  real(preal), dimension(3,3) :: vr, devnull, om_
601  integer :: ierr, i
602 
603  external :: dgeev
604 
605  om_ = om
606 
607  ! first get the rotation angle
608  t = 0.5_preal * (math_trace33(om) - 1.0_preal)
609  ax(4) = acos(math_clip(t,-1.0_preal,1.0_preal))
610 
611  if (deq0(ax(4))) then
612  ax(1:3) = [ 0.0_preal, 0.0_preal, 1.0_preal ]
613  else
614  call dgeev('N','V',3,om_,3,wr,wi,devnull,3,vr,3,work,size(work,1),ierr)
615  if (ierr /= 0) call io_error(401,ext_msg='Error in om2ax: DGEEV return not zero')
616 
617  i = maxloc(merge(1,0,ceq(cmplx(wr,wi,preal),cmplx(1.0_preal,0.0_preal,preal),tol=1.0e-14_preal)),dim=1)
618 
619 
620 
621  if (i == 0) call io_error(401,ext_msg='Error in om2ax Real: eigenvalue not found')
622  ax(1:3) = vr(1:3,i)
623  where ( dneq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
624  ax(1:3) = sign(ax(1:3),-p *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
625  endif
626 
627 end function om2ax
628 
629 
630 !---------------------------------------------------------------------------------------------------
633 !---------------------------------------------------------------------------------------------------
634 pure function om2ro(om) result(ro)
635 
636  real(preal), intent(in), dimension(3,3) :: om
637  real(preal), dimension(4) :: ro
638 
639  ro = eu2ro(om2eu(om))
640 
641 end function om2ro
642 
643 
644 !---------------------------------------------------------------------------------------------------
647 !---------------------------------------------------------------------------------------------------
648 function om2ho(om) result(ho)
649 
650  real(preal), intent(in), dimension(3,3) :: om
651  real(preal), dimension(3) :: ho
652 
653  ho = ax2ho(om2ax(om))
654 
655 end function om2ho
656 
657 
658 !---------------------------------------------------------------------------------------------------
661 !---------------------------------------------------------------------------------------------------
662 function om2cu(om) result(cu)
663 
664  real(preal), intent(in), dimension(3,3) :: om
665  real(preal), dimension(3) :: cu
666 
667  cu = ho2cu(om2ho(om))
668 
669 end function om2cu
670 
671 
672 !---------------------------------------------------------------------------------------------------
675 !---------------------------------------------------------------------------------------------------
676 pure function eu2qu(eu) result(qu)
677 
678  real(preal), intent(in), dimension(3) :: eu
679  real(preal), dimension(4) :: qu
680  real(preal), dimension(3) :: ee
681  real(preal) :: cphi, sphi
682 
683  ee = 0.5_preal*eu
684 
685  cphi = cos(ee(2))
686  sphi = sin(ee(2))
687 
688  qu = [ cphi*cos(ee(1)+ee(3)), &
689  -p*sphi*cos(ee(1)-ee(3)), &
690  -p*sphi*sin(ee(1)-ee(3)), &
691  -p*cphi*sin(ee(1)+ee(3))]
692  if(qu(1) < 0.0_preal) qu = qu * (-1.0_preal)
693 
694 end function eu2qu
695 
696 
697 !---------------------------------------------------------------------------------------------------
700 !---------------------------------------------------------------------------------------------------
701 pure function eu2om(eu) result(om)
702 
703  real(preal), intent(in), dimension(3) :: eu
704  real(preal), dimension(3,3) :: om
705 
706  real(preal), dimension(3) :: c, s
707 
708  c = cos(eu)
709  s = sin(eu)
710 
711  om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2)
712  om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2)
713  om(1,3) = s(3)*s(2)
714  om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2)
715  om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2)
716  om(2,3) = c(3)*s(2)
717  om(3,1) = s(1)*s(2)
718  om(3,2) = -c(1)*s(2)
719  om(3,3) = c(2)
720 
721  where(deq0(om)) om = 0.0_preal
722 
723 end function eu2om
724 
725 
726 !---------------------------------------------------------------------------------------------------
729 !---------------------------------------------------------------------------------------------------
730 pure function eu2ax(eu) result(ax)
731 
732  real(preal), intent(in), dimension(3) :: eu
733  real(preal), dimension(4) :: ax
734 
735  real(preal) :: t, delta, tau, alpha, sigma
736 
737  t = tan(eu(2)*0.5_preal)
738  sigma = 0.5_preal*(eu(1)+eu(3))
739  delta = 0.5_preal*(eu(1)-eu(3))
740  tau = sqrt(t**2+sin(sigma)**2)
741 
742  alpha = merge(pi, 2.0_preal*atan(tau/cos(sigma)), deq(sigma,pi*0.5_preal,tol=1.0e-15_preal))
743 
744  if (deq0(alpha)) then ! return a default identity axis-angle pair
745  ax = [ 0.0_preal, 0.0_preal, 1.0_preal, 0.0_preal ]
746  else
747  ax(1:3) = -p/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front
748  ax(4) = alpha
749  if (alpha < 0.0_preal) ax = -ax ! ensure alpha is positive
750  end if
751 
752 end function eu2ax
753 
754 
755 !---------------------------------------------------------------------------------------------------
758 !---------------------------------------------------------------------------------------------------
759 pure function eu2ro(eu) result(ro)
760 
761  real(preal), intent(in), dimension(3) :: eu
762  real(preal), dimension(4) :: ro
763 
764  ro = eu2ax(eu)
765  if (ro(4) >= pi) then
766  ro(4) = ieee_value(ro(4),ieee_positive_inf)
767  elseif(deq0(ro(4))) then
768  ro = [ 0.0_preal, 0.0_preal, p, 0.0_preal ]
769  else
770  ro(4) = tan(ro(4)*0.5_preal)
771  end if
772 
773 end function eu2ro
774 
775 
776 !---------------------------------------------------------------------------------------------------
779 !---------------------------------------------------------------------------------------------------
780 pure function eu2ho(eu) result(ho)
781 
782  real(preal), intent(in), dimension(3) :: eu
783  real(preal), dimension(3) :: ho
784 
785  ho = ax2ho(eu2ax(eu))
786 
787 end function eu2ho
788 
789 
790 !---------------------------------------------------------------------------------------------------
793 !---------------------------------------------------------------------------------------------------
794 function eu2cu(eu) result(cu)
795 
796  real(preal), intent(in), dimension(3) :: eu
797  real(preal), dimension(3) :: cu
798 
799  cu = ho2cu(eu2ho(eu))
800 
801 end function eu2cu
802 
803 
804 !---------------------------------------------------------------------------------------------------
807 !---------------------------------------------------------------------------------------------------
808 pure function ax2qu(ax) result(qu)
809 
810  real(preal), intent(in), dimension(4) :: ax
811  real(preal), dimension(4) :: qu
812 
813  real(preal) :: c, s
814 
815 
816  if (deq0(ax(4))) then
817  qu = [ 1.0_preal, 0.0_preal, 0.0_preal, 0.0_preal ]
818  else
819  c = cos(ax(4)*0.5_preal)
820  s = sin(ax(4)*0.5_preal)
821  qu = [ c, ax(1)*s, ax(2)*s, ax(3)*s ]
822  end if
823 
824 end function ax2qu
825 
826 
827 !---------------------------------------------------------------------------------------------------
830 !---------------------------------------------------------------------------------------------------
831 pure function ax2om(ax) result(om)
832 
833  real(preal), intent(in), dimension(4) :: ax
834  real(preal), dimension(3,3) :: om
835 
836  real(preal) :: q, c, s, omc
837 
838  c = cos(ax(4))
839  s = sin(ax(4))
840  omc = 1.0_preal-c
841 
842  om(1,1) = ax(1)**2*omc + c
843  om(2,2) = ax(2)**2*omc + c
844  om(3,3) = ax(3)**2*omc + c
845 
846  q = omc*ax(1)*ax(2)
847  om(1,2) = q + s*ax(3)
848  om(2,1) = q - s*ax(3)
849 
850  q = omc*ax(2)*ax(3)
851  om(2,3) = q + s*ax(1)
852  om(3,2) = q - s*ax(1)
853 
854  q = omc*ax(3)*ax(1)
855  om(3,1) = q + s*ax(2)
856  om(1,3) = q - s*ax(2)
857 
858  if (p > 0.0_preal) om = transpose(om)
859 
860 end function ax2om
861 
862 
863 !---------------------------------------------------------------------------------------------------
866 !---------------------------------------------------------------------------------------------------
867 pure function ax2eu(ax) result(eu)
868 
869  real(preal), intent(in), dimension(4) :: ax
870  real(preal), dimension(3) :: eu
871 
872  eu = om2eu(ax2om(ax))
873 
874 end function ax2eu
875 
876 
877 !---------------------------------------------------------------------------------------------------
880 !---------------------------------------------------------------------------------------------------
881 pure function ax2ro(ax) result(ro)
882 
883  real(preal), intent(in), dimension(4) :: ax
884  real(preal), dimension(4) :: ro
885 
886  real(preal), parameter :: thr = 1.0e-7_preal
887 
888  if (deq0(ax(4))) then
889  ro = [ 0.0_preal, 0.0_preal, p, 0.0_preal ]
890  else
891  ro(1:3) = ax(1:3)
892  ! we need to deal with the 180 degree case
893  ro(4) = merge(ieee_value(ro(4),ieee_positive_inf),tan(ax(4)*0.5_preal),abs(ax(4)-pi) < thr)
894  end if
895 
896 end function ax2ro
897 
898 
899 !---------------------------------------------------------------------------------------------------
902 !---------------------------------------------------------------------------------------------------
903 pure function ax2ho(ax) result(ho)
904 
905  real(preal), intent(in), dimension(4) :: ax
906  real(preal), dimension(3) :: ho
907 
908  real(preal) :: f
909 
910  f = 0.75_preal * ( ax(4) - sin(ax(4)) )
911  f = f**(1.0_preal/3.0_preal)
912  ho = ax(1:3) * f
913 
914 end function ax2ho
915 
916 
917 !---------------------------------------------------------------------------------------------------
920 !---------------------------------------------------------------------------------------------------
921 function ax2cu(ax) result(cu)
922 
923  real(preal), intent(in), dimension(4) :: ax
924  real(preal), dimension(3) :: cu
925 
926  cu = ho2cu(ax2ho(ax))
927 
928 end function ax2cu
929 
930 
931 !---------------------------------------------------------------------------------------------------
934 !---------------------------------------------------------------------------------------------------
935 pure function ro2qu(ro) result(qu)
936 
937  real(preal), intent(in), dimension(4) :: ro
938  real(preal), dimension(4) :: qu
939 
940  qu = ax2qu(ro2ax(ro))
941 
942 end function ro2qu
943 
944 
945 !---------------------------------------------------------------------------------------------------
948 !---------------------------------------------------------------------------------------------------
949 pure function ro2om(ro) result(om)
950 
951  real(preal), intent(in), dimension(4) :: ro
952  real(preal), dimension(3,3) :: om
953 
954  om = ax2om(ro2ax(ro))
955 
956 end function ro2om
957 
958 
959 !---------------------------------------------------------------------------------------------------
962 !---------------------------------------------------------------------------------------------------
963 pure function ro2eu(ro) result(eu)
964 
965  real(preal), intent(in), dimension(4) :: ro
966  real(preal), dimension(3) :: eu
967 
968  eu = om2eu(ro2om(ro))
969 
970 end function ro2eu
971 
972 
973 !---------------------------------------------------------------------------------------------------
976 !---------------------------------------------------------------------------------------------------
977 pure function ro2ax(ro) result(ax)
978 
979  real(preal), intent(in), dimension(4) :: ro
980  real(preal), dimension(4) :: ax
981 
982  real(preal) :: ta, angle
983 
984  ta = ro(4)
985 
986  if (.not. ieee_is_finite(ta)) then
987  ax = [ ro(1), ro(2), ro(3), pi ]
988  elseif (deq0(ta)) then
989  ax = [ 0.0_preal, 0.0_preal, 1.0_preal, 0.0_preal ]
990  else
991  angle = 2.0_preal*atan(ta)
992  ta = 1.0_preal/norm2(ro(1:3))
993  ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ]
994  end if
995 
996 end function ro2ax
997 
998 
999 !---------------------------------------------------------------------------------------------------
1002 !---------------------------------------------------------------------------------------------------
1003 pure function ro2ho(ro) result(ho)
1004 
1005  real(preal), intent(in), dimension(4) :: ro
1006  real(preal), dimension(3) :: ho
1007 
1008  real(preal) :: f
1009 
1010  if (deq0(norm2(ro(1:3)))) then
1011  ho = [ 0.0_preal, 0.0_preal, 0.0_preal ]
1012  else
1013  f = merge(2.0_preal*atan(ro(4)) - sin(2.0_preal*atan(ro(4))),pi, ieee_is_finite(ro(4)))
1014  ho = ro(1:3) * (0.75_preal*f)**(1.0_preal/3.0_preal)
1015  end if
1016 
1017 end function ro2ho
1018 
1019 
1020 !---------------------------------------------------------------------------------------------------
1023 !---------------------------------------------------------------------------------------------------
1024 pure function ro2cu(ro) result(cu)
1025 
1026  real(preal), intent(in), dimension(4) :: ro
1027  real(preal), dimension(3) :: cu
1028 
1029  cu = ho2cu(ro2ho(ro))
1030 
1031 end function ro2cu
1032 
1033 
1034 !---------------------------------------------------------------------------------------------------
1037 !---------------------------------------------------------------------------------------------------
1038 pure function ho2qu(ho) result(qu)
1039 
1040  real(preal), intent(in), dimension(3) :: ho
1041  real(preal), dimension(4) :: qu
1042 
1043  qu = ax2qu(ho2ax(ho))
1044 
1045 end function ho2qu
1046 
1047 
1048 !---------------------------------------------------------------------------------------------------
1051 !---------------------------------------------------------------------------------------------------
1052 pure function ho2om(ho) result(om)
1053 
1054  real(preal), intent(in), dimension(3) :: ho
1055  real(preal), dimension(3,3) :: om
1056 
1057  om = ax2om(ho2ax(ho))
1058 
1059 end function ho2om
1060 
1061 
1062 !---------------------------------------------------------------------------------------------------
1065 !---------------------------------------------------------------------------------------------------
1066 pure function ho2eu(ho) result(eu)
1067 
1068  real(preal), intent(in), dimension(3) :: ho
1069  real(preal), dimension(3) :: eu
1070 
1071  eu = ax2eu(ho2ax(ho))
1072 
1073 end function ho2eu
1074 
1075 
1076 !---------------------------------------------------------------------------------------------------
1079 !---------------------------------------------------------------------------------------------------
1080 pure function ho2ax(ho) result(ax)
1081 
1082  real(preal), intent(in), dimension(3) :: ho
1083  real(preal), dimension(4) :: ax
1084 
1085  integer :: i
1086  real(preal) :: hmag_squared, s, hm
1087  real(preal), parameter, dimension(16) :: &
1088  tfit = [ 1.0000000000018852_preal, -0.5000000002194847_preal, &
1089  -0.024999992127593126_preal, -0.003928701544781374_preal, &
1090  -0.0008152701535450438_preal, -0.0002009500426119712_preal, &
1091  -0.00002397986776071756_preal, -0.00008202868926605841_preal, &
1092  +0.00012448715042090092_preal, -0.0001749114214822577_preal, &
1093  +0.0001703481934140054_preal, -0.00012062065004116828_preal, &
1094  +0.000059719705868660826_preal, -0.00001980756723965647_preal, &
1095  +0.000003953714684212874_preal, -0.00000036555001439719544_preal ]
1096 
1097  ! normalize h and store the magnitude
1098  hmag_squared = sum(ho**2.0_preal)
1099  if (deq0(hmag_squared)) then
1100  ax = [ 0.0_preal, 0.0_preal, 1.0_preal, 0.0_preal ]
1101  else
1102  hm = hmag_squared
1103 
1104  ! convert the magnitude to the rotation angle
1105  s = tfit(1) + tfit(2) * hmag_squared
1106  do i=3,16
1107  hm = hm*hmag_squared
1108  s = s + tfit(i) * hm
1109  end do
1110  ax = [ho/sqrt(hmag_squared), 2.0_preal*acos(s)]
1111  end if
1112 
1113 end function ho2ax
1114 
1115 
1116 !---------------------------------------------------------------------------------------------------
1119 !---------------------------------------------------------------------------------------------------
1120 pure function ho2ro(ho) result(ro)
1121 
1122  real(preal), intent(in), dimension(3) :: ho
1123  real(preal), dimension(4) :: ro
1124 
1125  ro = ax2ro(ho2ax(ho))
1126 
1127 end function ho2ro
1128 
1129 
1130 !---------------------------------------------------------------------------------------------------
1133 !---------------------------------------------------------------------------------------------------
1134 pure function ho2cu(ho) result(cu)
1135 
1136  real(preal), intent(in), dimension(3) :: ho
1137  real(preal), dimension(3) :: cu
1138 
1139  cu = lambert_balltocube(ho)
1140 
1141 end function ho2cu
1142 
1143 
1144 !---------------------------------------------------------------------------------------------------
1147 !---------------------------------------------------------------------------------------------------
1148 pure function cu2qu(cu) result(qu)
1149 
1150  real(preal), intent(in), dimension(3) :: cu
1151  real(preal), dimension(4) :: qu
1152 
1153  qu = ho2qu(cu2ho(cu))
1154 
1155 end function cu2qu
1156 
1157 
1158 !---------------------------------------------------------------------------------------------------
1161 !---------------------------------------------------------------------------------------------------
1162 pure function cu2om(cu) result(om)
1163 
1164  real(preal), intent(in), dimension(3) :: cu
1165  real(preal), dimension(3,3) :: om
1166 
1167  om = ho2om(cu2ho(cu))
1168 
1169 end function cu2om
1170 
1171 
1172 !---------------------------------------------------------------------------------------------------
1175 !---------------------------------------------------------------------------------------------------
1176 pure function cu2eu(cu) result(eu)
1177 
1178  real(preal), intent(in), dimension(3) :: cu
1179  real(preal), dimension(3) :: eu
1180 
1181  eu = ho2eu(cu2ho(cu))
1182 
1183 end function cu2eu
1184 
1185 
1186 !---------------------------------------------------------------------------------------------------
1189 !---------------------------------------------------------------------------------------------------
1190 function cu2ax(cu) result(ax)
1191 
1192  real(preal), intent(in), dimension(3) :: cu
1193  real(preal), dimension(4) :: ax
1194 
1195  ax = ho2ax(cu2ho(cu))
1196 
1197 end function cu2ax
1198 
1199 
1200 !---------------------------------------------------------------------------------------------------
1203 !---------------------------------------------------------------------------------------------------
1204 pure function cu2ro(cu) result(ro)
1205 
1206  real(preal), intent(in), dimension(3) :: cu
1207  real(preal), dimension(4) :: ro
1208 
1209  ro = ho2ro(cu2ho(cu))
1210 
1211 end function cu2ro
1212 
1213 
1214 !---------------------------------------------------------------------------------------------------
1217 !---------------------------------------------------------------------------------------------------
1218 pure function cu2ho(cu) result(ho)
1219 
1220  real(preal), intent(in), dimension(3) :: cu
1221  real(preal), dimension(3) :: ho
1222 
1223  ho = lambert_cubetoball(cu)
1224 
1225 end function cu2ho
1226 
1227 
1228 !--------------------------------------------------------------------------------------------------
1230 !--------------------------------------------------------------------------------------------------
1231 subroutine unittest
1232 
1233  type(rotation) :: R
1234  real(pReal), dimension(4) :: qu, ax, ro
1235  real(pReal), dimension(3) :: x, eu, ho, v3
1236  real(pReal), dimension(3,3) :: om, t33
1237  real(pReal), dimension(3,3,3,3) :: t3333
1238  character(len=pStringLen) :: msg
1239  real :: A,B
1240  integer :: i
1241 
1242  do i=1,10
1243 
1244  msg = ''
1245 
1246 
1247  if(i<7) cycle
1248 
1249 
1250  if(i==1) then
1251  qu = om2qu(math_i3)
1252  elseif(i==2) then
1253  qu = eu2qu([0.0_preal,0.0_preal,0.0_preal])
1254  elseif(i==3) then
1255  qu = eu2qu([2.0_preal*pi,pi,2.0_preal*pi])
1256  elseif(i==4) then
1257  qu = [0.0_preal,0.0_preal,1.0_preal,0.0_preal]
1258  elseif(i==5) then
1259  qu = ro2qu([1.0_preal,0.0_preal,0.0_preal,ieee_value(1.0_preal, ieee_positive_inf)])
1260  elseif(i==6) then
1261  qu = ax2qu([1.0_preal,0.0_preal,0.0_preal,0.0_preal])
1262  else
1263  call random_number(x)
1264  a = sqrt(x(3))
1265  b = sqrt(1-0_preal -x(3))
1266  qu = [cos(2.0_preal*pi*x(1))*a,&
1267  sin(2.0_preal*pi*x(2))*b,&
1268  cos(2.0_preal*pi*x(2))*b,&
1269  sin(2.0_preal*pi*x(1))*a]
1270  if(qu(1)<0.0_preal) qu = qu * (-1.0_preal)
1271  endif
1272 
1273  if(dneq0(norm2(om2qu(qu2om(qu))-qu),1.0e-12_preal)) msg = trim(msg)//'om2qu/qu2om,'
1274  if(dneq0(norm2(eu2qu(qu2eu(qu))-qu),1.0e-12_preal)) msg = trim(msg)//'eu2qu/qu2eu,'
1275  if(dneq0(norm2(ax2qu(qu2ax(qu))-qu),1.0e-12_preal)) msg = trim(msg)//'ax2qu/qu2ax,'
1276  if(dneq0(norm2(ro2qu(qu2ro(qu))-qu),1.0e-12_preal)) msg = trim(msg)//'ro2qu/qu2ro,'
1277  if(dneq0(norm2(ho2qu(qu2ho(qu))-qu),1.0e-7_preal)) msg = trim(msg)//'ho2qu/qu2ho,'
1278  if(dneq0(norm2(cu2qu(qu2cu(qu))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2qu/qu2cu,'
1279 
1280 
1281  om = qu2om(qu)
1282 
1283  if(dneq0(norm2(om2qu(eu2om(om2eu(om)))-qu),1.0e-7_preal)) msg = trim(msg)//'eu2om/om2eu,'
1284  if(dneq0(norm2(om2qu(ax2om(om2ax(om)))-qu),1.0e-7_preal)) msg = trim(msg)//'ax2om/om2ax,'
1285  if(dneq0(norm2(om2qu(ro2om(om2ro(om)))-qu),1.0e-12_preal)) msg = trim(msg)//'ro2om/om2ro,'
1286  if(dneq0(norm2(om2qu(ho2om(om2ho(om)))-qu),1.0e-7_preal)) msg = trim(msg)//'ho2om/om2ho,'
1287  if(dneq0(norm2(om2qu(cu2om(om2cu(om)))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2om/om2cu,'
1288 
1289 
1290  eu = qu2eu(qu)
1291 
1292  if(dneq0(norm2(eu2qu(ax2eu(eu2ax(eu)))-qu),1.0e-12_preal)) msg = trim(msg)//'ax2eu/eu2ax,'
1293  if(dneq0(norm2(eu2qu(ro2eu(eu2ro(eu)))-qu),1.0e-12_preal)) msg = trim(msg)//'ro2eu/eu2ro,'
1294  if(dneq0(norm2(eu2qu(ho2eu(eu2ho(eu)))-qu),1.0e-7_preal)) msg = trim(msg)//'ho2eu/eu2ho,'
1295  if(dneq0(norm2(eu2qu(cu2eu(eu2cu(eu)))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2eu/eu2cu,'
1296 
1297 
1298  ax = qu2ax(qu)
1299 
1300  if(dneq0(norm2(ax2qu(ro2ax(ax2ro(ax)))-qu),1.0e-12_preal)) msg = trim(msg)//'ro2ax/ax2ro,'
1301  if(dneq0(norm2(ax2qu(ho2ax(ax2ho(ax)))-qu),1.0e-7_preal)) msg = trim(msg)//'ho2ax/ax2ho,'
1302  if(dneq0(norm2(ax2qu(cu2ax(ax2cu(ax)))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2ax/ax2cu,'
1303 
1304 
1305  ro = qu2ro(qu)
1306 
1307  if(dneq0(norm2(ro2qu(ho2ro(ro2ho(ro)))-qu),1.0e-7_preal)) msg = trim(msg)//'ho2ro/ro2ho,'
1308  if(dneq0(norm2(ro2qu(cu2ro(ro2cu(ro)))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2ro/ro2cu,'
1309 
1310 
1311  ho = qu2ho(qu)
1312 
1313  if(dneq0(norm2(ho2qu(cu2ho(ho2cu(ho)))-qu),1.0e-7_preal)) msg = trim(msg)//'cu2ho/ho2cu,'
1314 
1315 
1316  call r%fromMatrix(om)
1317 
1318  call random_number(v3)
1319  if(all(dneq(r%rotVector(r%rotVector(v3),active=.true.),v3,1.0e-12_preal))) &
1320  msg = trim(msg)//'rotVector,'
1321 
1322  call random_number(t33)
1323  if(all(dneq(r%rotTensor2(r%rotTensor2(t33),active=.true.),t33,1.0e-12_preal))) &
1324  msg = trim(msg)//'rotTensor2,'
1325 
1326  call random_number(t3333)
1327  if(all(dneq(r%rotTensor4(r%rotTensor4(t3333),active=.true.),t3333,1.0e-12_preal))) &
1328  msg = trim(msg)//'rotTensor4,'
1329 
1330  if(len_trim(msg) /= 0) call io_error(0,ext_msg=msg)
1331 
1332  enddo
1333 
1334 end subroutine unittest
1335 
1336 
1337 end module rotations
prec::deq
logical elemental pure function deq(a, b, tol)
equality comparison for float with double precision
Definition: prec.f90:123
quaternions::p
real(preal), parameter, public p
parameter for orientation conversion.
Definition: quaternions.f90:19
rotations
rotation storage and conversion
Definition: rotations.f90:53
lambert::lambert_balltocube
pure real(preal) function, dimension(3), public lambert_balltocube(xyz)
map from 3D ball to 3D cubic grid
Definition: Lambert.f90:132
rotations::ho2ax
pure real(preal) function, dimension(4) ho2ax(ho)
convert homochoric to axis angle pair
Definition: rotations.f90:1081
rotations::ax2eu
pure real(preal) function, dimension(3) ax2eu(ax)
convert axis angle pair to Euler angles
Definition: rotations.f90:868
rotations::om2ro
pure real(preal) function, dimension(4) om2ro(om)
convert rotation matrix to Rodrigues vector
Definition: rotations.f90:635
lambert
Mapping homochoric <-> cubochoric.
Definition: Lambert.f90:44
io::io_error
subroutine, public io_error(error_ID, el, ip, g, instance, ext_msg)
write error statements to standard out and terminate the Marc/spectral run with exit #9xxx
Definition: IO.f90:305
rotations::rotvector
pure real(preal) function, dimension(3) rotvector(self, v, active)
rotate a vector passively (default) or actively
Definition: rotations.f90:274
rotations::rottensor2
pure real(preal) function, dimension(3, 3) rottensor2(self, T, active)
rotate a rank-2 tensor passively (default) or actively
Definition: rotations.f90:311
rotations::rotation
Definition: rotations.f90:63
rotations::asaxisangle
pure real(preal) function, dimension(4) asaxisangle(self)
Definition: rotations.f90:128
rotations::asrodrigues
pure real(preal) function, dimension(4) asrodrigues(self)
Definition: rotations.f90:146
rotations::cu2om
pure real(preal) function, dimension(3, 3) cu2om(cu)
convert cubochoric to rotation matrix
Definition: rotations.f90:1163
rotations::fromaxisangle
subroutine fromaxisangle(self, ax, degrees, P)
Definition: rotations.f90:200
rotations::asmatrix
pure real(preal) function, dimension(3, 3) asmatrix(self)
Definition: rotations.f90:137
rotations::ho2qu
pure real(preal) function, dimension(4) ho2qu(ho)
convert homochoric to unit quaternion
Definition: rotations.f90:1039
math::math_sym3333to66
pure real(preal) function, dimension(6, 6) math_sym3333to66(m3333, weighted)
convert symmetric 3x3x3x3 matrix into 6x6 matrix
Definition: math.f90:778
rotations::fromquaternion
subroutine fromquaternion(self, qu)
Definition: rotations.f90:167
rotations::misorientation
pure elemental type(rotation) function misorientation(self, other)
misorientation
Definition: rotations.f90:391
rotations::cu2ro
pure real(preal) function, dimension(4) cu2ro(cu)
convert cubochoric to Rodrigues vector
Definition: rotations.f90:1205
math::inrad
real(preal), parameter inrad
conversion from degree into radian
Definition: math.f90:29
rotations::ax2ho
pure real(preal) function, dimension(3) ax2ho(ax)
convert axis angle pair to homochoric
Definition: rotations.f90:904
rotations::qu2eu
pure real(preal) function, dimension(3) qu2eu(qu)
convert unit quaternion to Euler angles
Definition: rotations.f90:435
rotations::frommatrix
subroutine frommatrix(self, om)
Definition: rotations.f90:230
rotations::qu2cu
pure real(preal) function, dimension(3) qu2cu(qu)
convert unit quaternion to cubochoric
Definition: rotations.f90:540
rotations::eu2ax
pure real(preal) function, dimension(4) eu2ax(eu)
convert euler to axis angle
Definition: rotations.f90:731
rotations::om2cu
real(preal) function, dimension(3) om2cu(om)
convert rotation matrix to cubochoric
Definition: rotations.f90:663
rotations::ashomochoric
pure real(preal) function, dimension(3) ashomochoric(self)
Definition: rotations.f90:155
rotations::rottensor4
pure real(preal) function, dimension(3, 3, 3, 3) rottensor4(self, T, active)
rotate a rank-4 tensor passively (default) or actively
Definition: rotations.f90:341
rotations::qu2om
pure real(preal) function, dimension(3, 3) qu2om(qu)
convert unit quaternion to rotation matrix
Definition: rotations.f90:405
rotations::qu2ho
pure real(preal) function, dimension(3) qu2ho(qu)
convert unit quaternion to homochoric
Definition: rotations.f90:516
rotations::fromeulers
subroutine fromeulers(self, eu, degrees)
Definition: rotations.f90:179
prec
setting precision for real and int type
Definition: prec.f90:13
quaternions::quaternion
Definition: quaternions.f90:21
math::math_det33
real(preal) pure function math_det33(m)
determinant of a 3x3 matrix
Definition: math.f90:623
quaternions::quaternions_init
subroutine, public quaternions_init
do self test
Definition: quaternions.f90:117
math::math_66tosym3333
pure real(preal) function, dimension(3, 3, 3, 3) math_66tosym3333(m66, weighted)
convert 66 matrix into symmetric 3x3x3x3 matrix
Definition: math.f90:806
rotations::cu2ax
real(preal) function, dimension(4) cu2ax(cu)
convert cubochoric to axis angle pair
Definition: rotations.f90:1191
rotations::rotrot__
pure elemental type(rotation) function rotrot__(self, R)
: Rotate a rotation
Definition: rotations.f90:247
rotations::qu2ro
pure real(preal) function, dimension(4) qu2ro(qu)
convert unit quaternion to Rodrigues vector
Definition: rotations.f90:489
rotations::ho2om
pure real(preal) function, dimension(3, 3) ho2om(ho)
convert homochoric to rotation matrix
Definition: rotations.f90:1053
rotations::rottensor4sym
pure real(preal) function, dimension(6, 6) rottensor4sym(self, T, active)
rotate a symmetric rank-4 tensor stored as (6,6) passively (default) or actively ToDo: Need to check ...
Definition: rotations.f90:372
quaternions::conjg
Definition: quaternions.f90:83
prec::dneq0
logical elemental pure function dneq0(a, tol)
inequality to 0 comparison for float with double precision
Definition: prec.f90:189
quaternions::abs
Definition: quaternions.f90:75
rotations::ax2om
pure real(preal) function, dimension(3, 3) ax2om(ax)
convert axis angle pair to orientation matrix
Definition: rotations.f90:832
rotations::om2eu
pure real(preal) function, dimension(3) om2eu(om)
orientation matrix to Euler angles
Definition: rotations.f90:569
io
input/output functions, partly depending on chosen solver
Definition: IO.f90:12
rotations::ax2cu
real(preal) function, dimension(3) ax2cu(ax)
convert axis angle pair to cubochoric
Definition: rotations.f90:922
rotations::ho2cu
pure real(preal) function, dimension(3) ho2cu(ho)
convert homochoric to cubochoric
Definition: rotations.f90:1135
rotations::ho2eu
pure real(preal) function, dimension(3) ho2eu(ho)
convert homochoric to Euler angles
Definition: rotations.f90:1067
rotations::om2qu
pure real(preal) function, dimension(4) om2qu(om)
convert rotation matrix to cubochoric
Definition: rotations.f90:555
prec::preal
integer, parameter preal
number with 15 significant digits, up to 1e+-307 (typically 64 bit)
Definition: prec.f90:20
rotations::eu2cu
real(preal) function, dimension(3) eu2cu(eu)
convert Euler angles to cubochoric
Definition: rotations.f90:795
prec::ceq
logical elemental pure function ceq(a, b, tol)
equality comparison for complex with double precision
Definition: prec.f90:210
math::math_clip
real(preal) pure elemental function math_clip(a, left, right)
limits a scalar value to a certain range (either one or two sided)
Definition: math.f90:1197
rotations::aseulers
pure real(preal) function, dimension(3) aseulers(self)
Definition: rotations.f90:119
prec::dneq
logical elemental pure function dneq(a, b, tol)
inequality comparison for float with double precision
Definition: prec.f90:146
rotations::om2ax
real(preal) function, dimension(4) om2ax(om)
convert orientation matrix to axis angle pair
Definition: rotations.f90:593
rotations::cu2eu
pure real(preal) function, dimension(3) cu2eu(cu)
convert cubochoric to Euler angles
Definition: rotations.f90:1177
quaternions
general quaternion math, not limited to unit quaternions
Definition: quaternions.f90:12
rotations::qu2ax
pure real(preal) function, dimension(4) qu2ax(qu)
convert unit quaternion to axis angle pair
Definition: rotations.f90:465
quaternions::real
Definition: quaternions.f90:95
math::math_trace33
real(preal) pure function math_trace33(m)
trace of a 3x3 matrix
Definition: math.f90:611
rotations::cu2qu
pure real(preal) function, dimension(4) cu2qu(cu)
convert cubochoric to unit quaternion
Definition: rotations.f90:1149
rotations::standardize
pure elemental subroutine standardize(self)
quaternion representation with positive q
Definition: rotations.f90:261
rotations::rotations_init
subroutine, public rotations_init
doing self test
Definition: rotations.f90:98
rotations::ro2om
pure real(preal) function, dimension(3, 3) ro2om(ro)
convert Rodrigues vector to rotation matrix
Definition: rotations.f90:950
rotations::ro2ax
pure real(preal) function, dimension(4) ro2ax(ro)
convert Rodrigues vector to axis angle pair
Definition: rotations.f90:978
rotations::ax2qu
pure real(preal) function, dimension(4) ax2qu(ax)
convert axis angle pair to quaternion
Definition: rotations.f90:809
rotations::asquaternion
pure real(preal) function, dimension(4) asquaternion(self)
Definition: rotations.f90:110
rotations::cu2ho
pure real(preal) function, dimension(3) cu2ho(cu)
convert cubochoric to homochoric
Definition: rotations.f90:1219
math
Mathematical library, including random number generation and tensor representations.
Definition: math.f90:12
rotations::ho2ro
pure real(preal) function, dimension(4) ho2ro(ho)
convert homochoric to Rodrigues vector
Definition: rotations.f90:1121
rotations::eu2ho
pure real(preal) function, dimension(3) eu2ho(eu)
convert Euler angles to homochoric
Definition: rotations.f90:781
rotations::eu2ro
pure real(preal) function, dimension(4) eu2ro(eu)
Euler angles to Rodrigues vector.
Definition: rotations.f90:760
rotations::om2ho
real(preal) function, dimension(3) om2ho(om)
convert rotation matrix to homochoric
Definition: rotations.f90:649
rotations::eu2om
pure real(preal) function, dimension(3, 3), public eu2om(eu)
Euler angles to orientation matrix.
Definition: rotations.f90:702
io::unittest
subroutine unittest
check correctness of some IO functions
Definition: IO.f90:662
rotations::ro2cu
pure real(preal) function, dimension(3) ro2cu(ro)
convert Rodrigues vector to cubochoric
Definition: rotations.f90:1025
rotations::ax2ro
pure real(preal) function, dimension(4) ax2ro(ax)
convert axis angle pair to Rodrigues vector
Definition: rotations.f90:882
lambert::lambert_cubetoball
pure real(preal) function, dimension(3), public lambert_cubetoball(cube)
map from 3D cubic grid to 3D ball
Definition: Lambert.f90:76
rotations::ro2qu
pure real(preal) function, dimension(4) ro2qu(ro)
convert Rodrigues vector to unit quaternion
Definition: rotations.f90:936
math::pi
real(preal), parameter pi
ratio of a circle's circumference to its diameter
Definition: math.f90:27
rotations::ro2ho
pure real(preal) function, dimension(3) ro2ho(ro)
convert Rodrigues vector to homochoric
Definition: rotations.f90:1004
rotations::ro2eu
pure real(preal) function, dimension(3) ro2eu(ro)
convert Rodrigues vector to Euler angles
Definition: rotations.f90:964
prec::deq0
logical elemental pure function deq0(a, tol)
equality to 0 comparison for float with double precision
Definition: prec.f90:166
math::math_i3
real(preal), dimension(3, 3), parameter math_i3
3x3 Identity
Definition: math.f90:32
rotations::eu2qu
pure real(preal) function, dimension(4) eu2qu(eu)
Euler angles to unit quaternion.
Definition: rotations.f90:677