grid_tools  1.4.0
 All Data Structures Files Functions Variables
pmat4.f90
Go to the documentation of this file.
1 
5 
45 module pmat4
46 !============================================================================
47 use pkind, only: spi,sp,dp,dpc
48 implicit none
49 private
50 public:: absv,normalized,orthogonalized, &
54  axtoq,qtoax, &
56  expmat,zntay,znfun, &
59 
60 interface absv; module procedure absv_s,absv_d; end interface
61 interface normalized;module procedure normalized_s,normalized_d;end interface
62 interface orthogonalized
63  module procedure orthogonalized_s,orthogonalized_d; end interface
64 interface cross_product
65  module procedure cross_product_s,cross_product_d, &
67 interface outer_product
69  end interface
70 interface triple_product
71  module procedure triple_product_s,triple_product_d; end interface
72 interface det; module procedure det_s,det_d,det_i,det_id; end interface
73 interface axial
74  module procedure axial3_s,axial3_d,axial33_s,axial33_d; end interface
75 interface diag
77  end interface
78 interface trace; module procedure trace_s,trace_d,trace_i; end interface
79 interface identity; module procedure identity_i,identity3_i; end interface
80 interface huarea; module procedure huarea_s,huarea_d; end interface
81 interface sarea
83  end interface
84 interface dlltoxy; module procedure dlltoxy_s,dlltoxy_d; end interface
85 interface hav; module procedure hav_s, hav_d; end interface
86 interface normalize;module procedure normalize_s,normalize_d; end interface
87 interface gram
89  end interface
90 interface rowops; module procedure rowops; end interface
91 interface corral; module procedure corral; end interface
92 interface rottoax; module procedure rottoax; end interface
93 interface axtorot; module procedure axtorot; end interface
94 interface spintoq; module procedure spintoq; end interface
95 interface qtospin; module procedure qtospin; end interface
96 interface rottoq; module procedure rottoq; end interface
97 interface qtorot; module procedure qtorot; end interface
98 interface axtoq; module procedure axtoq; end interface
99 interface qtoax; module procedure qtoax; end interface
100 interface setem; module procedure setem; end interface
101 interface mulqq; module procedure mulqq; end interface
102 interface expmat; module procedure expmat,expmatd,expmatdd; end interface
103 interface zntay; module procedure zntay; end interface
104 interface znfun; module procedure znfun; end interface
105 interface ctoz; module procedure ctoz; end interface
106 interface ztoc; module procedure ztoc,ztocd; end interface
107 interface setmobius;module procedure setmobius,zsetmobius; end interface
108 interface mobius; module procedure zmobius,cmobius; end interface
109 interface mobiusi; module procedure zmobiusi; end interface
110 
111 contains
112 
118 function absv_s(a)result(s)! [absv]
119 implicit none
120 real(sp),dimension(:),intent(in):: a
121 real(sp) :: s
122 s=sqrt(dot_product(a,a))
123 end function absv_s
124 
130 function absv_d(a)result(s)! [absv]
131 implicit none
132 real(dp),dimension(:),intent(in):: a
133 real(dp) :: s
134 s=sqrt(dot_product(a,a))
135 end function absv_d
136 
142 function normalized_s(a)result(b)! [normalized]
143 use pietc_s, only: u0
144 implicit none
145 real(sp),dimension(:),intent(IN):: a
146 real(sp),dimension(size(a)) :: b
147 real(sp) :: s
148 s=absv_s(a); if(s==u0)then; b=u0;else;b=a/s;endif
149 end function normalized_s
150 
156 function normalized_d(a)result(b)! [normalized]
157 use pietc, only: u0
158 implicit none
159 real(dp),dimension(:),intent(IN):: a
160 real(dp),dimension(size(a)) :: b
161 real(dp) :: s
162 s=absv_d(a); if(s==u0)then; b=u0;else;b=a/s;endif
163 end function normalized_d
164 
171 function orthogonalized_s(u,a)result(b)! [orthogonalized]
172 implicit none
173 real(sp),dimension(:),intent(in):: u,a
174 real(sp),dimension(size(u)) :: b
175 real(sp) :: s
176 ! Note: this routine assumes u is already normalized
177 s=dot_product(u,a); b=a-u*s
178 end function orthogonalized_s
179 
186 function orthogonalized_d(u,a)result(b)! [orthogonalized]
187 implicit none
188 real(dp),dimension(:),intent(in):: u,a
189 real(dp),dimension(size(u)) :: b
190 real(dp) :: s
191 ! Note: this routine assumes u is already normalized
192 s=dot_product(u,a); b=a-u*s
193 end function orthogonalized_d
194 
201 function cross_product_s(a,b)result(c)! [cross_product]
202 implicit none
203 real(sp),dimension(3),intent(in):: a,b
204 real(sp),dimension(3) :: c
205 c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1)
206 end function cross_product_s
207 
214 function cross_product_d(a,b)result(c)! [cross_product]
215 implicit none
216 real(dp),dimension(3),intent(in):: a,b
217 real(dp),dimension(3) :: c
218 c(1)=a(2)*b(3)-a(3)*b(2); c(2)=a(3)*b(1)-a(1)*b(3); c(3)=a(1)*b(2)-a(2)*b(1)
219 end function cross_product_d
220 
231 function triple_cross_product_s(u,v,w)result(x)! [cross_product]
232 implicit none
233 real(sp),dimension(4),intent(in ):: u,v,w
234 real(sp),dimension(4) :: x
235 real(sp):: uv12,uv13,uv14,uv23,uv24,uv34
236 uv12=u(1)*v(2)-u(2)*v(1); uv13=u(1)*v(3)-u(3)*v(1); uv14=u(1)*v(4)-u(4)*v(1)
237  uv23=u(2)*v(3)-u(3)*v(2); uv24=u(2)*v(4)-u(4)*v(2)
238  uv34=u(3)*v(4)-u(4)*v(3)
239 x(1)=-uv23*w(4)+uv24*w(3)-uv34*w(2)
240 x(2)= uv13*w(4)-uv14*w(3) +uv34*w(1)
241 x(3)=-uv12*w(4) +uv14*w(2)-uv24*w(1)
242 x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1)
243 end function triple_cross_product_s
244 
252 function triple_cross_product_d(u,v,w)result(x)! [cross_product]
253 implicit none
254 real(dp),dimension(4),intent(in ):: u,v,w
255 real(dp),dimension(4) :: x
256 real(dp):: uv12,uv13,uv14,uv23,uv24,uv34
257 uv12=u(1)*v(2)-u(2)*v(1); uv13=u(1)*v(3)-u(3)*v(1); uv14=u(1)*v(4)-u(4)*v(1)
258  uv23=u(2)*v(3)-u(3)*v(2); uv24=u(2)*v(4)-u(4)*v(2)
259  uv34=u(3)*v(4)-u(4)*v(3)
260 x(1)=-uv23*w(4)+uv24*w(3)-uv34*w(2)
261 x(2)= uv13*w(4)-uv14*w(3) +uv34*w(1)
262 x(3)=-uv12*w(4) +uv14*w(2)-uv24*w(1)
263 x(4)= uv12*w(3)-uv13*w(2)+uv23*w(1)
264 end function triple_cross_product_d
265 
272 function outer_product_s(a,b)result(c)! [outer_product]
273 implicit none
274 real(sp),dimension(:), intent(in ):: a
275 real(sp),dimension(:), intent(in ):: b
276 real(sp),DIMENSION(size(a),size(b)):: c
277 integer(spi) :: nb,i
278 nb=size(b)
279 do i=1,nb; c(:,i)=a*b(i); enddo
280 end function outer_product_s
281 
288 function outer_product_d(a,b)result(c)! [outer_product]
289 implicit none
290 real(dp),dimension(:), intent(in ):: a
291 real(dp),dimension(:), intent(in ):: b
292 real(dp),dimension(size(a),size(b)):: c
293 integer(spi) :: nb,i
294 nb=size(b)
295 do i=1,nb; c(:,i)=a*b(i); enddo
296 end function outer_product_d
297 
304 function outer_product_i(a,b)result(c)! [outer_product]
305 implicit none
306 integer(spi),dimension(:), intent(in ):: a
307 integer(spi),dimension(:), intent(in ):: b
308 integer(spi),dimension(size(a),size(b)):: c
309 integer(spi) :: nb,i
310 nb=size(b)
311 do i=1,nb; c(:,i)=a*b(i); enddo
312 end function outer_product_i
313 
321 function triple_product_s(a,b,c)result(tripleproduct)! [triple_product]
322 implicit none
323 real(sp),dimension(3),intent(IN ):: a,b,c
324 real(sp) :: tripleproduct
325 tripleproduct=dot_product( cross_product(a,b),c )
326 end function triple_product_s
327 
335 function triple_product_d(a,b,c)result(tripleproduct)! [triple_product]
336 implicit none
337 real(dp),dimension(3),intent(IN ):: a,b,c
338 real(dp) :: tripleproduct
339 tripleproduct=dot_product( cross_product(a,b),c )
340 end function triple_product_d
341 
347 function det_s(a)result(det)! [det]
348 use pietc_s, only: u0
349 implicit none
350 real(sp),dimension(:,:),intent(IN ) :: a
351 real(sp) :: det
352 real(sp),dimension(size(a,1),size(a,1)):: b
353 integer(spi) :: n,nrank
354 n=size(a,1)
355 if(n==3)then
356  det=triple_product(a(:,1),a(:,2),a(:,3))
357 else
358  call gram(a,b,nrank,det)
359  if(nrank<n)det=u0
360 endif
361 end function det_s
362 
368 function det_d(a)result(det)! [det]
369 use pietc, only: u0
370 implicit none
371 real(dp),dimension(:,:),intent(IN ) :: a
372 real(dp) :: det
373 real(dp),dimension(size(a,1),size(a,1)):: b
374 integer(spi) :: n,nrank
375 n=size(a,1)
376 if(n==3)then
377  det=triple_product(a(:,1),a(:,2),a(:,3))
378 else
379  call gram(a,b,nrank,det)
380  if(nrank<n)det=u0
381 endif
382 end function det_d
383 
389 function det_i(a)result(idet)! [det]
390 implicit none
391 integer(spi), dimension(:,:),intent(IN ):: a
392 integer(spi) :: idet
393 real(dp),dimension(size(a,1),size(a,2)):: b
394 real(dp) :: bdet
395 b=a; bdet=det(b); idet=nint(bdet)
396 end function det_i
397 
403 function det_id(a)result(idet)! [det]
404 use pkind, only: dp,dpi
405 implicit none
406 integer(dpi), dimension(:,:),intent(IN ):: a
407 integer(dpi) :: idet
408 real(dp),dimension(size(a,1),size(a,2)) :: b
409 real(dp) :: bdet
410 b=a; bdet=det(b); idet=nint(bdet)
411 end function det_id
412 
419 function axial3_s(a)result(b)! [axial]
420 use pietc_s, only: u0
421 implicit none
422 real(sp),dimension(3),intent(IN ):: a
423 real(sp),dimension(3,3) :: b
424 b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3)
425 end function axial3_s
426 
433 function axial3_d(a)result(b)! [axial]
434 use pietc, only: u0
435 implicit none
436 real(dp),dimension(3),intent(IN ):: a
437 real(dp),dimension(3,3) :: b
438 b=u0;b(3,2)=a(1);b(1,3)=a(2);b(2,1)=a(3);b(2,3)=-a(1);b(3,1)=-a(2);b(1,2)=-a(3)
439 end function axial3_d
440 
447 function axial33_s(b)result(a)! [axial]
448 use pietc_s, only: o2
449 implicit none
450 real(sp),dimension(3,3),intent(IN ):: b
451 real(sp),dimension(3) :: a
452 a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2
453 end function axial33_s
454 
461 function axial33_d(b)result(a)! [axial]
462 use pietc, only: o2
463 implicit none
464 real(dp),dimension(3,3),intent(IN ):: b
465 real(dp),dimension(3) :: a
466 a(1)=(b(3,2)-b(2,3))*o2; a(2)=(b(1,3)-b(3,1))*o2; a(3)=(b(2,1)-b(1,2))*o2
467 end function axial33_d
468 
475 function diagn_s(a)result(b)! [diag]
476 use pietc, only: u0
477 implicit none
478 real(sp),dimension(:),intent(IN ) :: a
479 real(sp),dimension(size(a),size(a)):: b
480 integer(spi) :: n,i
481 n=size(a)
482 b=u0; do i=1,n; b(i,i)=a(i); enddo
483 end function diagn_s
484 
491 function diagn_d(a)result(b)! [diag]
492 use pietc, only: u0
493 implicit none
494 real(dp),dimension(:),intent(IN ) :: a
495 real(dp),dimension(size(a),size(a)):: b
496 integer(spi) :: n,i
497 n=size(a)
498 b=u0; do i=1,n; b(i,i)=a(i); enddo
499 end function diagn_d
500 
507 function diagn_i(a)result(b)! [diag]
508 implicit none
509 integer(spi),dimension(:),intent(IN ) :: a
510 integer(spi),dimension(size(a),size(a)):: b
511 integer(spi) :: n,i
512 n=size(a)
513 b=0; do i=1,n; b(i,i)=a(i); enddo
514 end function diagn_i
515 
522 function diagnn_s(b)result(a)! [diag]
523 implicit none
524 real(sp),dimension(:,:),intent(IN ):: b
525 real(sp),dimension(size(b,1)) :: a
526 integer(spi) :: n,i
527 n=size(b,1)
528 do i=1,n; a(i)=b(i,i); enddo
529 end function diagnn_s
530 
537 function diagnn_d(b)result(a)! [diag]
538 implicit none
539 real(dp),dimension(:,:),intent(IN ):: b
540 real(dp),dimension(size(b,1)) :: a
541 integer(spi) :: n,i
542 n=size(b,1)
543 do i=1,n; a(i)=b(i,i); enddo
544 end function diagnn_d
545 
552 function diagnn_i(b)result(a)! [diag]
553 implicit none
554 integer(spi),dimension(:,:),intent(IN ):: b
555 integer(spi),dimension(size(b,1)) :: a
556 integer(spi) :: n,i
557 n=size(b,1)
558 do i=1,n; a(i)=b(i,i); enddo
559 end function diagnn_i
560 
566 function trace_s(b)result(s)! [trace]
567 implicit none
568 real(sp),dimension(:,:),intent(IN ):: b
569 real(sp) :: s
570 s=sum(diag(b))
571 end function trace_s
572 
578 function trace_d(b)result(s)! [trace]
579 implicit none
580 real(dp),dimension(:,:),intent(IN ):: b
581 real(dp) :: s
582 s=sum(diag(b))
583 end function trace_d
584 
590 function trace_i(b)result(s)! [trace]
591 implicit none
592 integer(spi),dimension(:,:),intent(IN ):: b
593 integer(spi) :: s
594 s=sum(diag(b))
595 end function trace_i
596 
602 function identity_i(n)result(a)! [identity]
603 implicit none
604 integer(spi),intent(IN ) :: n
605 integer(spi),dimension(n,n):: a
606 integer(spi) :: i
607 a=0; do i=1,n; a(i,i)=1; enddo
608 end function identity_i
609 
614 function identity3_i()result(a)! [identity]
615 implicit none
616 integer(spi),dimension(3,3):: a
617 integer(spi) :: i
618 a=0; do i=1,3; a(i,i)=1; enddo
619 end function identity3_i
620 
629 function huarea_s(sa,sb)result(area)! [huarea]
630 implicit none
631 real(sp),intent(IN ):: sa,sb
632 real(sp) :: area
633 real(sp) :: ca,cb
634 ca=sqrt(1-sa**2)
635 cb=sqrt(1-sb**2)
636 area=asin(sa*sb/(1+ca*cb))
637 end function huarea_s
638 
647 function huarea_d(sa,sb)result(area)! [huarea]
648 implicit none
649 real(dp),intent(IN ):: sa,sb
650 real(dp) :: area
651 real(dp) :: ca,cb
652 ca=sqrt(1-sa**2)
653 cb=sqrt(1-sb**2)
654 area=asin(sa*sb/(1+ca*cb))
655 end function huarea_d
656 
667 function sarea_s(v1,v2,v3)result(area)! [sarea]
668 use pietc_s, only: zero=>u0
669 implicit none
670 real(sp),dimension(3),intent(IN ):: v1,v2,v3
671 real(sp) :: area
672 real(sp) :: s123,a1,a2,b,d1,d2,d3
673 real(sp),dimension(3) :: u0,u1,u2,u3,x,y
674 area=zero
675 u1=normalized(v1); u2=normalized(v2); u3=normalized(v3)
676 s123=triple_product(u1,u2,u3)
677 if(s123==zero)return
678 
679 d1=dot_product(u3-u2,u3-u2)
680 d2=dot_product(u1-u3,u1-u3)
681 d3=dot_product(u2-u1,u2-u1)
682 
683 ! Triangle that is not degenerate. Cyclically permute, so side 3 is longest:
684 if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
685 if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
686 y=normalized( cross_product(u1,u2) )
687 b=dot_product(y,u3)
688 u0=normalized( u3-y*b )
689 x=cross_product(y,u0)
690 a1=-dot_product(x,u1-u0); a2= dot_product(x,u2-u0)
691 area=huarea(a1,b)+huarea(a2,b)
692 
693 contains
694 
704 subroutine cyclic(u1,u2,u3,d1,d2,d3)
705 implicit none
706  real(sp),dimension(3),intent(INOUT):: u1,u2,u3
707  real(sp), intent(INOUT):: d1,d2,d3
708  real(sp),dimension(3) :: ut
709  real(sp) :: dt
710  dt=d1; d1=d2; d2=d3; d3=dt
711  ut=u1; u1=u2; u2=u3; u3=ut
712  end subroutine cyclic
713 end function sarea_s
714 
722 function sarea_d(v1,v2,v3)result(area)! [sarea]
723 use pietc, only: zero=>u0
724 implicit none
725 real(dp),dimension(3),intent(IN ):: v1,v2,v3
726 real(dp) :: area
727 real(dp) :: s123,a1,a2,b,d1,d2,d3
728 real(dp),dimension(3) :: u0,u1,u2,u3,x,y
729 area=zero
730 u1=normalized(v1); u2=normalized(v2); u3=normalized(v3)
731 s123=triple_product(u1,u2,u3)
732 if(s123==zero)return
733 
734 d1=dot_product(u3-u2,u3-u2)
735 d2=dot_product(u1-u3,u1-u3)
736 d3=dot_product(u2-u1,u2-u1)
737 
738 ! Triangle that is not degenerate. Cyclically permute, so side 3 is longest:
739 if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
740 if(d3<d1 .or. d3<d2)call cyclic(u1,u2,u3,d1,d2,d3)
741 y=normalized( cross_product(u1,u2) )
742 b=dot_product(y,u3)
743 u0=normalized( u3-y*b )
744 x=cross_product(y,u0)
745 a1=-dot_product(x,u1-u0); a2= dot_product(x,u2-u0)
746 area=huarea(a1,b)+huarea(a2,b)
747 
748 contains
749 
759 subroutine cyclic(u1,u2,u3,d1,d2,d3)
760 implicit none
761  real(dp),dimension(3),intent(INOUT):: u1,u2,u3
762  real(dp), intent(INOUT):: d1,d2,d3
763  real(dp),dimension(3) :: ut
764  real(dp) :: dt
765  dt=d1; d1=d2; d2=d3; d3=dt
766  ut=u1; u1=u2; u2=u3; u3=ut
767  end subroutine cyclic
768 end function sarea_d
769 
785 function dtarea_s(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea]
786 use pietc_s, only: u0,u1
787 implicit none
788 real(sp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb
789 real(sp) :: area
790 real(sp),dimension(2):: x2a,x2b,xb,yb
791 real(sp) :: sb,ssb,cb,xa,sa,ca,sc,cc
792 call dlltoxy(rlat,drlata,drlona,x2a)
793 call dlltoxy(rlat,drlatb,drlonb,x2b)
794 ssb=dot_product(x2b,x2b); sb=sqrt(ssb)
795 if(sb==u0)then; area=u0; return; endif
796 cb=sqrt(u1-ssb)
797 ! Construct 2D normalized right-handed basis vectors with xb pointing to B:
798 xb=x2b/sb
799 yb=(/-xb(2),xb(1)/)
800 xa=dot_product(xb,x2a)
801 sa=dot_product(yb,x2a)
802 ca=sqrt(u1-sa**2)
803 sc=xa/ca
804 cc=sqrt(u1-sc**2)
805 sb=sb*cc-cb*sc
806 area=huarea(-sa,sb)+huarea(sc,-sa)
807 end function dtarea_s
808 
824 function dtarea_d(rlat,drlata,drlona,drlatb,drlonb) result(area)! [sarea]
825 use pietc, only: u0,u1
826 implicit none
827 real(dp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb
828 real(dp) :: area
829 real(dp),dimension(2):: x2a,x2b,xb,yb
830 real(dp) :: sb,ssb,cb,xa,sa,ca,sc,cc
831 call dlltoxy(rlat,drlata,drlona,x2a)
832 call dlltoxy(rlat,drlatb,drlonb,x2b)
833 ssb=dot_product(x2b,x2b); sb=sqrt(ssb)
834 if(sb==u0)then; area=u0; return; endif
835 cb=sqrt(u1-ssb)
836 ! Construct 2D normalized right-handed basis vectors with xb pointing to B:
837 xb=x2b/sb
838 yb=(/-xb(2),xb(1)/)
839 xa=dot_product(xb,x2a)
840 sa=dot_product(yb,x2a)
841 ca=sqrt(u1-sa**2)
842 sc=xa/ca
843 cc=sqrt(u1-sc**2)
844 sb=sb*cc-cb*sc
845 area=huarea(-sa,sb)+huarea(sc,-sa)
846 end function dtarea_d
847 
867 function dqarea_s &! [sarea]
868  (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area)
869 implicit none
870 real(sp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc
871 real(sp) :: area
872 area=sarea(rlat,drlata,drlona,drlatb,drlonb)&
873  -sarea(rlat,drlatc,drlonc,drlatb,drlonb)
874 end function dqarea_s
875 
876 
896 function dqarea_d &! [sarea]
897  (rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc) result(area)
898 implicit none
899 real(dp),intent(in ):: rlat,drlata,drlona,drlatb,drlonb,drlatc,drlonc
900 real(dp) :: area
901 area=sarea(rlat,drlata,drlona,drlatb,drlonb)&
902  -sarea(rlat,drlatc,drlonc,drlatb,drlonb)
903 end function dqarea_d
904 
915 subroutine dlltoxy_s(rlat,drlat,drlon,x2)! [dlltoxy]
916 use pietc_s, only: u2
917 implicit none
918 real(sp), intent(in ):: rlat,drlat,drlon
919 real(sp),dimension(2),intent(out):: x2
920 real(sp):: clata
921 clata=cos(rlat+drlat)
922 x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/)
923 end subroutine dlltoxy_s
924 
935 subroutine dlltoxy_d(rlat,drlat,drlon,x2)! [dlltoxy]
936 use pietc, only: u2
937 implicit none
938 real(dp), intent(in ):: rlat,drlat,drlon
939 real(dp),dimension(2),intent(out):: x2
940 real(dp):: clata
941 clata=cos(rlat+drlat)
942 x2=(/clata*sin(drlon),sin(drlat)+u2*sin(rlat)*clata*hav(drlon)/)
943 end subroutine dlltoxy_d
944 
950 function hav_s(t) result(a)! [hav]
951 use pietc_s, only: o2
952 implicit none
953 real(sp),intent(in ):: t
954 real(sp) :: a
955 a=(sin(t*o2))**2
956 end function hav_s
957 
963 function hav_d(t) result(a)! [hav]
964 use pietc, only: o2
965 implicit none
966 real(dp),intent(in ):: t
967 real(dp) :: a
968 a=(sin(t*o2))**2
969 end function hav_d
970 
975 subroutine normalize_s(v)! [normalize]
976 use pietc_s, only: u0,u1
977 implicit none
978 real(sp),dimension(:),intent(inout):: v
979 real(sp) :: s
980 s=absv(v); if(s==0)then; v=u0; v(1)=u1; else; v=v/s; endif
981 end subroutine normalize_s
982 
987 subroutine normalize_d(v)! [normalize]
988 use pietc, only: u0,u1
989 implicit none
990 real(dp),dimension(:),intent(inout):: v
991 real(dp) :: s
992 s=absv(v); if(s==u0)then; v=0; v(1)=u1; else; v=v/s; endif
993 end subroutine normalize_d
994 
1006 subroutine gram_s(as,b,nrank,det)! [gram]
1007 use pietc_s, only: u0,u1
1008 implicit none
1009 real(sp),dimension(:,:),intent(IN ) :: as
1010 real(sp),dimension(:,:),intent(OUT) :: b
1011 integer(spi), intent(OUT) :: nrank
1012 real(sp), intent(OUT) :: det
1013 real(sp),parameter :: crit=1.e-5_sp
1014 real(sp),dimension(size(as,1),size(as,2)):: a
1015 real(sp),dimension(size(as,2),size(as,1)):: ab
1016 real(sp),dimension(size(as,1)) :: tv,w
1017 real(sp) :: val,s,vcrit
1018 integer(spi) :: i,j,k,l,m,n
1019 integer(spi),dimension(2) :: ii
1020 n=size(as,1)
1021 m=size(as,2)
1022 if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions'
1023 a=as
1024 b=identity(n)
1025 det=u1
1026 val=maxval(abs(a))
1027 if(val==u0)then
1028  nrank=0
1029  return
1030 endif
1031 vcrit=val*crit
1032 nrank=min(n,m)
1033 do k=1,n
1034  if(k>nrank)exit
1035  ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) )
1036  ii =maxloc( abs( ab(k:m,k:n)) )+k-1
1037  val=maxval( abs( ab(k:m,k:n)) )
1038  if(val<=vcrit)then
1039  nrank=k-1
1040  exit
1041  endif
1042  i=ii(1)
1043  j=ii(2)
1044  tv=b(:,j)
1045  b(:,j)=-b(:,k)
1046  b(:,k)=tv
1047  tv=a(:,i)
1048  a(:,i)=-a(:,k)
1049  a(:,k)=tv
1050  w(k:n)=matmul( transpose(b(:,k:n)),tv )
1051  b(:,k)=matmul(b(:,k:n),w(k:n) )
1052  s=dot_product(b(:,k),b(:,k))
1053  s=sqrt(s)
1054  if(w(k)<u0)s=-s
1055  det=det*s
1056  b(:,k)=b(:,k)/s
1057  do l=k,n
1058  do j=l+1,n
1059  s=dot_product(b(:,l),b(:,j))
1060  b(:,j)=normalized( b(:,j)-b(:,l)*s )
1061  enddo
1062  enddo
1063 enddo
1064 end subroutine gram_s
1065 
1077 subroutine gram_d(as,b,nrank,det)! [gram]
1078 use pietc, only: u0,u1
1079 implicit none
1080 real(dp),dimension(:,:),intent(IN ) :: as
1081 real(dp),dimension(:,:),intent(OUT) :: b
1082 integer(spi), intent(OUT) :: nrank
1083 real(dp), intent(OUT) :: det
1084 real(dp),parameter :: crit=1.e-9_dp
1085 real(dp),dimension(size(as,1),size(as,2)):: a
1086 real(dp),dimension(size(as,2),size(as,1)):: ab
1087 real(dp),dimension(size(as,1)) :: tv,w
1088 real(dp) :: val,s,vcrit
1089 integer(spi) :: i,j,k,l,m,n
1090 integer(spi),dimension(2) :: ii
1091 n=size(as,1)
1092 m=size(as,2)
1093 if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions'
1094 a=as
1095 b=identity(n)
1096 det=u1
1097 val=maxval(abs(a))
1098 if(val==u0)then
1099  nrank=0
1100  return
1101 endif
1102 vcrit=val*crit
1103 nrank=min(n,m)
1104 do k=1,n
1105  if(k>nrank)exit
1106  ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) )
1107  ii =maxloc( abs( ab(k:m,k:n)) )+k-1
1108  val=maxval( abs( ab(k:m,k:n)) )
1109  if(val<=vcrit)then
1110  nrank=k-1
1111  exit
1112  endif
1113  i=ii(1)
1114  j=ii(2)
1115  tv=b(:,j)
1116  b(:,j)=-b(:,k)
1117  b(:,k)=tv
1118  tv=a(:,i)
1119  a(:,i)=-a(:,k)
1120  a(:,k)=tv
1121  w(k:n)=matmul( transpose(b(:,k:n)),tv )
1122  b(:,k)=matmul(b(:,k:n),w(k:n) )
1123  s=dot_product(b(:,k),b(:,k))
1124  s=sqrt(s)
1125  if(w(k)<u0)s=-s
1126  det=det*s
1127  b(:,k)=b(:,k)/s
1128  do l=k,n
1129  do j=l+1,n
1130  s=dot_product(b(:,l),b(:,j))
1131  b(:,j)=normalized( b(:,j)-b(:,l)*s )
1132  enddo
1133  enddo
1134 enddo
1135 end subroutine gram_d
1136 
1149 subroutine graml_d(as,b,nrank,detsign,ldet)! [gram]
1150 use pietc, only: u0
1151 implicit none
1152 real(dp),dimension(:,:),intent(IN ) :: as
1153 real(dp),dimension(:,:),intent(OUT) :: b
1154 integer(spi), intent(OUT) :: nrank
1155 integer(spi), intent(out) :: detsign
1156 real(dp), intent(OUT) :: ldet
1157 real(dp),parameter :: crit=1.e-9_dp
1158 real(dp),dimension(size(as,1),size(as,2)):: a
1159 real(dp),dimension(size(as,2),size(as,1)):: ab
1160 real(dp),dimension(size(as,1)) :: tv,w
1161 real(dp) :: val,s,vcrit
1162 integer(spi) :: i,j,k,l,m,n
1163 integer(spi),dimension(2) :: ii
1164 detsign=1
1165 n=size(as,1)
1166 m=size(as,2)
1167 if(n/=size(b,1) .or. n/=size(b,2))stop 'In gram; incompatible dimensions'
1168 a=as
1169 b=identity(n)
1170 ldet=u0
1171 val=maxval(abs(a))
1172 if(val==u0)then
1173  nrank=0
1174  return
1175 endif
1176 vcrit=val*crit
1177 nrank=min(n,m)
1178 do k=1,n
1179  if(k>nrank)exit
1180  ab(k:m,k:n)=matmul( transpose(a(:,k:m)),b(:,k:n) )
1181  ii =maxloc( abs( ab(k:m,k:n)) )+k-1
1182  val=maxval( abs( ab(k:m,k:n)) )
1183  if(val<=vcrit)then
1184  nrank=k-1
1185  exit
1186  endif
1187  i=ii(1)
1188  j=ii(2)
1189  tv=b(:,j)
1190  b(:,j)=-b(:,k)
1191  b(:,k)=tv
1192  tv=a(:,i)
1193  a(:,i)=-a(:,k)
1194  a(:,k)=tv
1195  w(k:n)=matmul( transpose(b(:,k:n)),tv )
1196  b(:,k)=matmul(b(:,k:n),w(k:n) )
1197  s=dot_product(b(:,k),b(:,k))
1198  s=sqrt(s)
1199  if(w(k)<u0)s=-s
1200  if(s<0)then
1201  ldet=ldet+log(-s)
1202  detsign=-detsign
1203  elseif(s>u0)then
1204  ldet=ldet+log(s)
1205  else
1206  detsign=0
1207  endif
1208 
1209  b(:,k)=b(:,k)/s
1210  do l=k,n
1211  do j=l+1,n
1212  s=dot_product(b(:,l),b(:,j))
1213  b(:,j)=normalized( b(:,j)-b(:,l)*s )
1214  enddo
1215  enddo
1216 enddo
1217 end subroutine graml_d
1218 
1219 
1226 subroutine plaingram_s(b,nrank)! [gram]
1227 use pietc_s, only: u0
1228 implicit none
1229 real(sp),dimension(:,:),intent(INOUT) :: b
1230 integer(spi), intent( OUT) :: nrank
1231 real(sp),parameter :: crit=1.e-5_sp
1232 real(sp) :: val,vcrit
1233 integer(spi) :: j,k,n
1234 n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square'
1235 val=maxval(abs(b))
1236 nrank=0
1237 if(val==0)then
1238  b=u0
1239  return
1240 endif
1241 vcrit=val*crit
1242 do k=1,n
1243  val=sqrt(dot_product(b(:,k),b(:,k)))
1244  if(val<=vcrit)then
1245  b(:,k:n)=u0
1246  return
1247  endif
1248  b(:,k)=b(:,k)/val
1249  nrank=k
1250  do j=k+1,n
1251  b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j))
1252  enddo
1253 enddo
1254 end subroutine plaingram_s
1255 
1262 subroutine plaingram_d(b,nrank)! [gram]
1263 use pietc, only: u0
1264 implicit none
1265 real(dp),dimension(:,:),intent(INOUT):: b
1266 integer(spi), intent( OUT):: nrank
1267 real(dp),parameter:: crit=1.e-9_dp
1268 real(dp) :: val,vcrit
1269 integer(spi) :: j,k,n
1270 n=size(b,1); if(n/=size(b,2))stop 'In gram; matrix needs to be square'
1271 val=maxval(abs(b))
1272 nrank=0
1273 if(val==u0)then
1274  b=u0
1275  return
1276 endif
1277 vcrit=val*crit
1278 do k=1,n
1279  val=sqrt(dot_product(b(:,k),b(:,k)))
1280  if(val<=vcrit)then
1281  b(:,k:n)=u0
1282  return
1283  endif
1284  b(:,k)=b(:,k)/val
1285  nrank=k
1286  do j=k+1,n
1287  b(:,j)=b(:,j)-b(:,k)*dot_product(b(:,k),b(:,j))
1288  enddo
1289 enddo
1290 end subroutine plaingram_d
1291 
1314 subroutine rowgram(m,n,a,ipiv,tt,b,rank)! [gram]
1315 use pietc, only: u0,u1
1316 implicit none
1317 integer(spi), intent(IN ):: m,n
1318 real(dp),dimension(m,n), intent(in ):: a
1319 integer(spi),dimension(n),intent(out):: ipiv
1320 real(dp),dimension(m,n), intent(out):: tt
1321 real(dp),dimension(n,n), intent(out):: b
1322 integer(spi), intent(out):: rank
1323 real(dp),parameter :: eps=1.e-13_dp,epss=eps**2
1324 real(dp),dimension(m,n) :: aa
1325 real(dp),dimension(n) :: rowv
1326 real(dp),dimension(m) :: p
1327 real(dp) :: maxp,nepss
1328 integer(spi),dimension(1):: jloc
1329 integer(spi) :: i,ii,iii,j,maxi
1330 if(m<n)stop 'In rowgram; this routines needs m>=n please'
1331 nepss=n*epss
1332 rank=n
1333 aa=a
1334 tt=u0
1335 do ii=1,n
1336 ! At this stage, all rows less than ii are already orthonormalized and are
1337 ! orthogonal to all rows at and beyond ii. Find the norms of these lower
1338 ! rows and pivot the largest of them into position ii:
1339  maxp=u0
1340  maxi=ii
1341  do i=ii,m
1342  p(i)=dot_product(aa(i,:),aa(i,:))
1343  if(p(i)>maxp)then
1344  maxp=p(i)
1345  maxi=i
1346  endif
1347  enddo
1348  if(maxp<nepss)then
1349  b=u0
1350  b(1:ii-1,:)=aa(1:ii-1,:)
1351 ! fill the remaining rows, ii:n, of b with remaining orthonormal rows at random
1352  do iii=ii,n
1353 ! find the column of b for which the maximum element is the smallest:
1354  do j=1,n
1355  rowv(j)=maxval(abs(b(1:iii-1,j)))
1356  enddo
1357  jloc=minloc(rowv)
1358  j=jloc(1)
1359  b(iii,j)=u1
1360  do i=1,iii-1
1361  maxp=dot_product(b(i,:),b(iii,:))
1362  b(iii,:)=b(iii,:)-b(i,:)*maxp
1363  enddo
1364  maxp=sqrt(dot_product(b(iii,:),b(iii,:)))
1365  b(iii,:)=b(iii,:)/maxp
1366  enddo
1367  rank=ii-1
1368  return
1369  endif
1370  ipiv(ii)=maxi
1371  if(maxi/=ii)then
1372  rowv =aa(ii, :)
1373  aa(ii, :)=aa(maxi,:)
1374  aa(maxi,:)=rowv
1375  endif
1376  maxp=sqrt(maxp)
1377  tt(ii,ii)=maxp
1378  aa(ii,:)=aa(ii,:)/maxp
1379 ! Adjust all rows below to make them orthogonal to new row ii
1380  do i=ii+1,m
1381  maxp=dot_product(aa(ii,:),aa(i,:))
1382  tt(i,ii)=maxp
1383  aa(i,:)=aa(i,:)-aa(ii,:)*maxp
1384  enddo
1385 enddo
1386 b=aa(1:n,:)
1387 end subroutine rowgram
1388 
1399 subroutine rowops(m,n,ipiv,tt,v,vv)! [rowops]
1400 implicit none
1401 integer(spi), intent(in ):: m,n
1402 integer(spi),dimension(n),intent(in ):: ipiv
1403 real(dp),dimension(m,n), intent(in ):: tt
1404 real(dp),dimension(m), intent(in ):: v
1405 real(dp),dimension(m), intent(out):: vv
1406 integer(spi):: i,j,k
1407 real(dp) :: p
1408 vv=v
1409 do j=1,n
1410  k=ipiv(j)
1411  if(k/=j)then
1412  p=vv(j)
1413  vv(j)=vv(k)
1414  vv(k)=p
1415  endif
1416  vv(j)=vv(j)/tt(j,j)
1417  do i=j+1,m
1418  vv(i)=vv(i)-vv(j)*tt(i,j)
1419  enddo
1420 enddo
1421 end subroutine rowops
1422 
1440 subroutine corral(m,n,mask,a,d,aa,e)! [corral]
1441 use pietc, only: u0,u1
1442 use pmat, only: inv
1443 implicit none
1444 integer(spi), intent(in ):: m,n
1445 logical, dimension(m,n),intent(in ):: mask
1446 real(dp),dimension(m,n),intent(in ):: a
1447 real(dp),dimension(m ),intent(out):: d
1448 real(dp),dimension(m,n),intent(out):: aa
1449 real(dp),dimension( n),intent(out):: e
1450 real(dp),dimension(0:m+n,0:m+n):: g
1451 real(dp),dimension(0:m+n) :: h
1452 integer(spi) :: i,j,k,nh
1453 nh=1+m+n
1454 aa=u0
1455 do j=1,n
1456 do i=1,m
1457  if(mask(i,j))aa(i,j)=log(abs(a(i,j)))
1458 enddo
1459 enddo
1460 h=u0
1461 g=u0
1462 ! Equations on row 0 enforcing Lagrange multiplier F-constraint:
1463 do j=1,n
1464  k=m+j
1465  g(0,k)=u1
1466 enddo
1467 ! Equations on rows 1:m minimizing row sums of quadratic terms:
1468 do i=1,m
1469  do j=1,n
1470  k=m+j
1471  if(mask(i,j))then
1472  g(i,i)=g(i,i)-u1
1473  g(i,k)=-u1
1474  h(i)=h(i)-aa(i,j)
1475  endif
1476  enddo
1477 enddo
1478 ! Equations on rows m+1:m+n minimizing col sums subject to constraint
1479 do j=1,n
1480  k=m+j
1481  g(k,0)=u1
1482  do i=1,m
1483  if(mask(i,j))then
1484  g(k,k)=g(k,k)-u1
1485  g(k,i)=-u1
1486  h(k)=h(k)-aa(i,j)
1487  endif
1488  enddo
1489 enddo
1490 ! Invert the normal equations:
1491 call inv(g,h)
1492 ! Exponentiate the parts that become final scaling diagnonal matrices d and e:
1493 do i=1,m
1494  d(i)=exp(h(i))
1495 enddo
1496 do j=1,n
1497  k=m+j
1498  e(j)=exp(h(k))
1499 enddo
1500 ! Compute the rescaled matrix directly:
1501 do j=1,n
1502 do i=1,m
1503  aa(i,j)=a(i,j)/(d(i)*e(j))
1504 enddo
1505 enddo
1506 end subroutine corral
1507 
1517 subroutine rottoax(orth33,ax3)! [rottoax]
1518 implicit none
1519 real(dp),dimension(3,3),intent(in ):: orth33
1520 real(dp),dimension(3), intent(out):: ax3
1521 real(dp),dimension(3,3) :: plane
1522 real(dp),dimension(3) :: x,y,z
1523 real(dp) :: s
1524 integer(spi),dimension(1):: ii
1525 integer(spi) :: i,j,k
1526 plane=orth33-identity()! Columns must be coplanar vectors
1527 do i=1,3; z(i)=dot_product(plane(:,i),plane(:,i)); enddo
1528 ii=minloc(z)
1529 k=ii(1); i=1+mod(k,3); j=1+mod(i,3)
1530 ax3=cross_product(plane(:,i),plane(:,j))
1531 s=absv(ax3); if(s==0)return
1532 ax3=ax3/s ! <- temporarily a unit vector pointing along rotation axis
1533 ! Construct a unit 2D basis, x,y, in the plane of rotation
1534 x=normalized(cross_product(ax3,plane(:,j)))
1535 y=cross_product(ax3,x)
1536 z=matmul(orth33,x)! Rotate x by the original rotation matrix
1537 ax3=ax3*atan2(dot_product(y,z),dot_product(x,z))! multiply ax3 by the angle
1538 end subroutine rottoax
1539 
1547 subroutine axtorot(ax3,orth33)! [axtorot]
1548 implicit none
1549 real(dp),dimension(3), intent(in ):: ax3
1550 real(dp),dimension(3,3),intent(out):: orth33
1551 real(dp),dimension(3,3):: ax33
1552 real(dp) :: d
1553 ax33=axial(ax3); call expmat(3,ax33,orth33,d)
1554 end subroutine axtorot
1555 
1561 subroutine spintoq(cspin,q)! [spintoq]
1562 implicit none
1563 complex(dpc),dimension(2,2),intent(IN ):: cspin
1564 real(dp), dimension(0:3),intent(OUT):: q
1565 q(0)=real(cspin(1,1)); q(3)=aimag(cspin(1,1))
1566 q(2)=real(cspin(2,1)); q(1)=aimag(cspin(2,1))
1567 end subroutine spintoq
1568 
1574 subroutine qtospin(q,cspin)! [qtospin]
1575 implicit none
1576 real(dp), dimension(0:3),intent(IN ):: q
1577 complex(dpc),dimension(2,2),intent(OUT):: cspin
1578 cspin(1,1)=cmplx( q(0), q(3))
1579 cspin(2,1)=cmplx( q(2), q(1))
1580 cspin(1,2)=cmplx(-q(2), q(1))
1581 cspin(2,2)=cmplx( q(0),-q(3))
1582 end subroutine qtospin
1583 
1589 subroutine rottoq(rot,q)! [rottoq]
1590 use pietc, only: zero=>u0,one=>u1,two=>u2
1591 implicit none
1592 real(dp),dimension(3,3),intent(IN ):: rot
1593 real(dp),dimension(0:3),intent(OUT):: q
1594 real(dp),dimension(3,3) :: t1,t2
1595 real(dp),dimension(3) :: u1,u2
1596 real(dp) :: gamma,gammah,s,ss
1597 integer(spi) :: i,j
1598 integer(spi),dimension(1):: ii
1599 ! construct the orthogonal matrix, t1, whose third row is the rotation axis
1600 ! of rot:
1601 t1=rot; do i=1,3; t1(i,i)=t1(i,i)-1; u1(i)=dot_product(t1(i,:),t1(i,:)); enddo
1602 ii=maxloc(u1); j=ii(1); ss=u1(j)
1603 if(ss<1.e-16_dp)then
1604  q=zero; q(0)=one; return
1605 endif
1606 t1(j,:)=t1(j,:)/sqrt(ss)
1607 if(j/=1)then
1608  u2 =t1(1,:)
1609  t1(1,:)=t1(j,:)
1610  t1(j,:)=u2
1611 endif
1612 do i=2,3
1613  t1(i,:)=t1(i,:)-dot_product(t1(1,:),t1(i,:))*t1(1,:)
1614  u1(i)=dot_product(t1(i,:),t1(i,:))
1615 enddo
1616 if(u1(3)>u1(2))then
1617  j=3
1618 else
1619  j=2
1620 endif
1621 ss=u1(j)
1622 if(ss==zero)stop 'In rotov; invalid rot'
1623 if(j/=2)t1(2,:)=t1(3,:)
1624 t1(2,:)=t1(2,:)/sqrt(ss)
1625 ! Form t1(3,:) as the cross product of t1(1,:) and t1(2,:)
1626 t1(3,1)=t1(1,2)*t1(2,3)-t1(1,3)*t1(2,2)
1627 t1(3,2)=t1(1,3)*t1(2,1)-t1(1,1)*t1(2,3)
1628 t1(3,3)=t1(1,1)*t1(2,2)-t1(1,2)*t1(2,1)
1629 ! Project rot into the frame whose axes are the rows of t1:
1630 t2=matmul(t1,matmul(rot,transpose(t1)))
1631 ! Obtain the rotation angle, gamma, implied by rot, and gammah=gamma/2:
1632 gamma=atan2(t2(2,1),t2(1,1)); gammah=gamma/two
1633 ! Hence deduce coefficients (in the form of a real 4-vector) of one of the two
1634 ! possible equivalent spinors:
1635 s=sin(gammah)
1636 q(0)=cos(gammah)
1637 q(1:3)=t1(3,:)*s
1638 end subroutine rottoq
1639 
1645 subroutine qtorot(q,rot)! [qtorot]
1646 implicit none
1647 real(dp),dimension(0:3),intent(IN ):: q
1648 real(dp),dimension(3,3),intent(OUT):: rot
1649 call setem(q(0),q(1),q(2),q(3),rot)
1650 end subroutine qtorot
1651 
1657 subroutine axtoq(v,q)! [axtoq]
1658 implicit none
1659 real(dp),dimension(3), intent(in ):: v
1660 real(dp),dimension(0:3),intent(out):: q
1661 real(dp),dimension(3,3):: rot
1662 call axtorot(v,rot)
1663 call rottoq(rot,q)
1664 end subroutine axtoq
1665 
1671 subroutine qtoax(q,v)! [qtoax]
1672 implicit none
1673 real(dp),dimension(0:3),intent(in ):: q
1674 real(dp),dimension(3), intent(out):: v
1675 real(dp),dimension(3,3):: rot
1676 call qtorot(q,rot)
1677 call rottoax(rot,v)
1678 end subroutine qtoax
1679 
1689 subroutine setem(c,d,e,g,r)! [setem]
1690 implicit none
1691 real(dp), intent(IN ):: c,d,e,g
1692 real(dp),dimension(3,3),intent(OUT):: r
1693 real(dp):: cc,dd,ee,gg,de,dg,eg,dc,ec,gc
1694 cc=c*c; dd=d*d; ee=e*e; gg=g*g
1695 de=d*e; dg=d*g; eg=e*g
1696 dc=d*c; ec=e*c; gc=g*c
1697 r(1,1)=cc+dd-ee-gg; r(2,2)=cc-dd+ee-gg; r(3,3)=cc-dd-ee+gg
1698 r(2,3)=2*(eg-dc); r(3,1)=2*(dg-ec); r(1,2)=2*(de-gc)
1699 r(3,2)=2*(eg+dc); r(1,3)=2*(dg+ec); r(2,1)=2*(de+gc)
1700 end subroutine setem
1701 
1708 function mulqq(a,b)result(c)! [mulqq]
1709 implicit none
1710 real(dp),dimension(0:3),intent(IN ):: a,b
1711 real(dp),dimension(0:3) :: c
1712 c(0)=a(0)*b(0) -a(1)*b(1) -a(2)*b(2) -a(3)*b(3)
1713 c(1)=a(0)*b(1) +a(1)*b(0) +a(2)*b(3) -a(3)*b(2)
1714 c(2)=a(0)*b(2) +a(2)*b(0) +a(3)*b(1) -a(1)*b(3)
1715 c(3)=a(0)*b(3) +a(3)*b(0) +a(1)*b(2) -a(2)*b(1)
1716 end function mulqq
1717 
1728 subroutine expmat(n,a,b,detb)! [expmat]
1729 use pietc, only: u0,u1,u2,o2
1730 implicit none
1731 integer(spi), intent(IN ):: n
1732 real(dp),dimension(n,n),intent(IN ):: a
1733 real(dp),dimension(n,n),intent(OUT):: b
1734 real(dp), intent(OUT):: detb
1735 integer(spi),parameter :: l=5
1736 real(dp),dimension(n,n):: c,p
1737 real(dp) :: t
1738 integer(spi) :: i,m
1739 m=10+floor(log(u1+maxval(abs(a)))/log(u2))
1740 t=o2**m
1741 c=a*t
1742 p=c
1743 b=p
1744 do i=2,l
1745  p=matmul(p,c)/i
1746  b=b+p
1747 enddo
1748 do i=1,m
1749  b=b*u2+matmul(b,b)
1750 enddo
1751 do i=1,n
1752  b(i,i)=b(i,i)+u1
1753 enddo
1754 detb=u0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb)
1755 end subroutine expmat
1756 
1766 subroutine expmatd(n,a,b,bd,detb,detbd)! [expmat]
1767 use pietc, only: u0,u1,u2,o2
1768 implicit none
1769 integer(spi), intent(IN ):: n
1770 real(dp),dimension(n,n), intent(IN ):: a
1771 real(dp),dimension(n,n), intent(OUT):: b
1772 real(dp),dimension(n,n,(n*(n+1))/2),intent(OUT):: bd
1773 real(dp), intent(OUT):: detb
1774 real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd
1775 integer(spi),parameter :: l=5
1776 real(dp),dimension(n,n) :: c,p
1777 real(dp),dimension(n,n,(n*(n+1))/2):: pd,cd
1778 real(dp) :: t
1779 integer(spi) :: i,j,k,m,n1
1780 n1=(n*(n+1))*o2
1781 m=10+floor(log(u1+maxval(abs(a)))/log(u2))
1782 t=o2**m
1783 c=a*t
1784 p=c
1785 pd=u0
1786 do k=1,n
1787  pd(k,k,k)=t
1788 enddo
1789 k=n
1790 do i=1,n-1
1791  do j=i+1,n
1792  k=k+1
1793  pd(i,j,k)=t
1794  pd(j,i,k)=t
1795  enddo
1796 enddo
1797 if(k/=n1)stop 'In expmatd; n1 is inconsistent with n'
1798 cd=pd
1799 b=p
1800 bd=pd
1801 do i=2,l
1802  do k=1,n1
1803  pd(:,:,k)=(matmul(cd(:,:,k),p)+matmul(c,pd(:,:,k)))/i
1804  enddo
1805  p=matmul(c,p)/i
1806  b=b+p
1807  bd=bd+pd
1808 enddo
1809 do i=1,m
1810  do k=1,n1
1811  bd(:,:,k)=2*bd(:,:,k)+matmul(bd(:,:,k),b)+matmul(b,bd(:,:,k))
1812  enddo
1813  b=b*u2+matmul(b,b)
1814 enddo
1815 do i=1,n
1816  b(i,i)=b(i,i)+u1
1817 enddo
1818 detb=u0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb)
1819 detbd=u0; do k=1,n; detbd(k)=detb; enddo
1820 end subroutine expmatd
1821 
1833 subroutine expmatdd(n,a,b,bd,bdd,detb,detbd,detbdd)! [expmat]
1834 use pietc, only: u0,u1,u2,o2
1835 implicit none
1836 integer(spi), intent(IN ):: n
1837 real(dp),dimension(n,n), intent(IN ):: a
1838 real(dp),dimension(n,n), intent(OUT):: b
1839 real(dp),dimension(n,n,(n*(n+1))/2), intent(OUT):: bd
1840 real(dp),dimension(n,n,(n*(n+1))/2,(n*(n+1))/2),intent(OUT):: bdd
1841 real(dp), intent(OUT):: detb
1842 real(dp),dimension((n*(n+1))/2), intent(OUT):: detbd
1843 real(dp),dimension((n*(n+1))/2,(n*(n+1))/2), intent(OUT):: detbdd
1844 integer(spi),parameter :: l=5
1845 real(dp),dimension(n,n) :: c,p
1846 real(dp),dimension(n,n,(n*(n+1))/2) :: pd,cd
1847 real(dp),dimension(n,n,(n*(n+1))/2,(n*(n+1))/2):: pdd,cdd
1848 real(dp) :: t
1849 integer(spi) :: i,j,k,ki,kj,m,n1
1850 n1=(n*(n+1))/2
1851 m=10+floor(log(u1+maxval(abs(a)))/log(u2))
1852 t=o2**m
1853 c=a*t
1854 p=c
1855 pd=u0
1856 pdd=u0
1857 do k=1,n
1858  pd(k,k,k)=t
1859 enddo
1860 k=n
1861 do i=1,n-1
1862  do j=i+1,n
1863  k=k+1
1864  pd(i,j,k)=t
1865  pd(j,i,k)=t
1866  enddo
1867 enddo
1868 if(k/=n1)stop 'In expmatd; n1 is inconsistent with n'
1869 cd=pd
1870 cdd=u0
1871 b=p
1872 bd=pd
1873 bdd=u0
1874 do i=2,l
1875  do ki=1,n1
1876  do kj=1,n1
1877  pdd(:,:,ki,kj)=(matmul(cd(:,:,ki),pd(:,:,kj)) &
1878  + matmul(cd(:,:,kj),pd(:,:,ki)) &
1879  + matmul(c,pdd(:,:,ki,kj)))/i
1880  enddo
1881  enddo
1882  do k=1,n1
1883  pd(:,:,k)=(matmul(cd(:,:,k),p)+matmul(c,pd(:,:,k)))/i
1884  enddo
1885  p=matmul(c,p)/i
1886  b=b+p
1887  bd=bd+pd
1888  bdd=bdd+pdd
1889 enddo
1890 do i=1,m
1891  do ki=1,n1
1892  do kj=1,n1
1893  bdd(:,:,ki,kj)=u2*bdd(:,:,ki,kj) &
1894  +matmul(bdd(:,:,ki,kj),b) &
1895  +matmul(bd(:,:,ki),bd(:,:,kj)) &
1896  +matmul(bd(:,:,kj),bd(:,:,ki)) &
1897  +matmul(b,bdd(:,:,ki,kj))
1898  enddo
1899  enddo
1900  do k=1,n1
1901  bd(:,:,k)=2*bd(:,:,k)+matmul(bd(:,:,k),b)+matmul(b,bd(:,:,k))
1902  enddo
1903  b=b*u2+matmul(b,b)
1904 enddo
1905 do i=1,n
1906  b(i,i)=b(i,i)+u1
1907 enddo
1908 detb=u0; do i=1,n; detb=detb+a(i,i); enddo; detb=exp(detb)
1909 detbd=u0; do k=1,n; detbd(k)=detb; enddo
1910 detbdd=u0; do ki=1,n; do kj=1,n; detbdd(ki,kj)=detb; enddo; enddo
1911 end subroutine expmatdd
1912 
1920 subroutine zntay(n,z,zn)! [zntay]
1921 use pietc, only: u2
1922 implicit none
1923 integer(spi), intent(IN ):: n
1924 real(dp), intent(IN ):: z
1925 real(dp), intent(OUT):: zn
1926 integer(spi),parameter:: ni=100
1927 real(dp),parameter :: eps0=1.e-16_dp
1928 integer(spi) :: i,i2,n2
1929 real(dp) :: t,eps,z2
1930 z2=z*u2
1931 n2=n*2
1932 t=1
1933 do i=1,n
1934  t=t/(i*2-1)
1935 enddo
1936 eps=t*eps0
1937 zn=t
1938 t=t
1939 do i=1,ni
1940  i2=i*2
1941  t=t*z2/(i2*(i2+n2-1))
1942  zn=zn+t
1943  if(abs(t)<eps)return
1944 enddo
1945 print'("In zntay; full complement of iterations used")'
1946 end subroutine zntay
1947 
1959 subroutine znfun(n,z,zn,znd,zndd,znddd)! [znfun]
1960 use pietc, only: u0,u2,u3
1961 implicit none
1962 integer(spi),intent(IN ):: n
1963 real(dp), intent(IN ):: z
1964 real(dp), intent(OUT):: zn,znd,zndd,znddd
1965 real(dp) :: z2,rz2
1966 integer(spi):: i,i2p3
1967 z2=abs(z*u2)
1968 rz2=sqrt(z2)
1969 if(z2<u2)then
1970  call zntay(n ,z,zn)
1971  call zntay(n+1,z,znd)
1972  call zntay(n+2,z,zndd)
1973  call zntay(n+3,z,znddd)
1974 else
1975  if(z>u0)then
1976  zn=cosh(rz2)
1977  znd=sinh(rz2)/rz2
1978  zndd=(zn-znd)/z2
1979  znddd=(znd-u3*zndd)/z2
1980  do i=1,n
1981  i2p3=i*2+3
1982  zn=znd
1983  znd=zndd
1984  zndd=znddd
1985  znddd=(znd-i2p3*zndd)/z2
1986  enddo
1987  else
1988  zn=cos(rz2)
1989  znd=sin(rz2)/rz2
1990  zndd=-(zn-znd)/z2
1991  znddd=-(znd-u3*zndd)/z2
1992  do i=1,n
1993  i2p3=i*2+3
1994  zn=znd
1995  znd=zndd
1996  zndd=znddd
1997  znddd=-(znd-i2p3*zndd)/z2
1998  enddo
1999  endif
2000 endif
2001 end subroutine znfun
2002 
2018 
2026 subroutine ctoz(v, z,infz)! [ctoz]
2027 use pietc, only: u0,u1
2028 implicit none
2029 real(dp),dimension(3),intent(IN ):: v
2030 complex(dpc), intent(OUT):: z
2031 logical, intent(OUT):: infz
2032 real(dp) :: rr,zzpi
2033 infz=.false.
2034 z=cmplx(v(1),v(2),dpc)
2035 if(v(3)>u0)then
2036  zzpi=u1/(u1+v(3))
2037 else
2038  rr=v(1)**2+v(2)**2
2039  infz=(rr==u0); if(infz)return ! <- The point is mapped to infinity (90S)
2040  zzpi=(u1-v(3))/rr
2041 endif
2042 z=z*zzpi
2043 end subroutine ctoz
2044 
2052 subroutine ztoc(z,infz, v)! [ztoc]
2053 implicit none
2054 complex(dpc), intent(IN ):: z
2055 logical, intent(IN ):: infz
2056 real(dp),dimension(3),intent(OUT):: v
2057 real(dp),parameter:: zero=0_dp,one=1_dp,two=2_dp
2058 real(dp) :: r,q,rs,rsc,rsbi
2059 if(infz)then; v=(/zero,zero,-one/); return; endif
2060 r=real(z); q=aimag(z); rs=r*r+q*q
2061 rsc=one-rs
2062 rsbi=one/(one+rs)
2063 v(1)=two*rsbi*r
2064 v(2)=two*rsbi*q
2065 v(3)=rsc*rsbi
2066 end subroutine ztoc
2067 
2082 subroutine ztocd(z,infz, v,vd)! [ztoc]
2083 implicit none
2084 complex(dpc), intent(IN ):: z
2085 logical, intent(IN ):: infz
2086 real(dp),dimension(3), intent(OUT):: v
2087 complex(dpc),dimension(3),intent(OUT):: vd
2088 real(dp),parameter :: zero=0_dp,one=1_dp,two=2_dp,four=4_dp
2089 real(dp) :: r,q,rs,rsc,rsbi,rsbis
2090 real(dp),dimension(3):: u1,u2
2091 integer(spi) :: i
2092 if(infz)then; v=(/zero,zero,-one/); return; endif
2093 r=real(z); q=aimag(z); rs=r*r+q*q
2094 rsc=one-rs
2095 rsbi=one/(one+rs)
2096 rsbis=rsbi**2
2097 v(1)=two*rsbi*r
2098 v(2)=two*rsbi*q
2099 v(3)=rsc*rsbi
2100 u1(1)=two*(one+q*q-r*r)*rsbis
2101 u1(2)=-four*r*q*rsbis
2102 u1(3)=-four*r*rsbis
2103 u2=cross_product(v,u1)
2104 do i=1,3
2105  vd(i)=cmplx(u1(i),-u2(i),dpc)
2106 enddo
2107 end subroutine ztocd
2108 
2122 subroutine setmobius(xc0,xc1,xc2, aa,bb,cc,dd)! [setmobius]
2123 implicit none
2124 real(dp),dimension(3),intent(IN ):: xc0,xc1,xc2
2125 complex(dpc), intent(OUT):: aa,bb,cc,dd
2126 real(dp),parameter:: zero=0_dp,one=1_dp
2127 logical :: infz0,infz1,infz2
2128 complex(dpc) :: z0,z1,z2,z02,z10,z21
2129 call ctoz(xc0,z0,infz0)
2130 call ctoz(xc1,z1,infz1)
2131 call ctoz(xc2,z2,infz2)
2132 z21=z2-z1
2133 z02=z0-z2
2134 z10=z1-z0
2135 
2136 if( (z0==z1.and.infz0.eqv.infz1).or.&
2137  (z1==z2.and.infz1.eqv.infz2).or.&
2138  (z2==z0.and.infz2.eqv.infz0)) &
2139  stop 'In setmobius; anchor points must be distinct'
2140 if(infz2 .or. (.not.infz0 .and. abs(z0)<abs(z2)))then
2141 ! z0 is finite and smaller than z2:
2142  if(infz1)then
2143  aa=one/sqrt(z02) ! <- z1 is infinite
2144  elseif(infz2)then
2145  aa=one/sqrt(z10) ! <- z2 is infinite
2146  else
2147  aa=sqrt(-z21/(z02*z10)) ! <- all zs are finite
2148  endif
2149  bb=-z0*aa
2150  if(infz1)then
2151  cc=aa ! <- z1 is infinite
2152  dd=-z2*aa !
2153  elseif(infz2)then
2154  cc=zero ! <- z2 is infinite
2155  dd=z10*aa !
2156  else
2157  cc=-(z10/z21)*aa ! <- all zs are finite
2158  dd= z2*(z10/z21)*aa !
2159  endif
2160 else
2161 ! z2 is finite and smaller than z0:
2162  if(infz0)then
2163  cc=one/sqrt(z21) ! <- z0 is inifinite
2164  elseif(infz1)then
2165  cc=one/sqrt(z02) ! <- z1 is infinite
2166  else
2167  cc=sqrt(-z10/(z02*z21)) ! <- all zs are finite
2168  endif
2169  dd=-z2*cc
2170  if(infz0)then
2171  aa=zero ! <- z0 is inifinite
2172  bb=-z21*cc !
2173  elseif(infz1)then
2174  aa=cc ! <- z1 is infinite
2175  bb=-z0*cc !
2176  else
2177  aa=(-z21/z10)*cc ! <- all zs are finite
2178  bb=z0*(z21/z10)*cc !
2179  endif
2180 endif
2181 end subroutine setmobius
2182 
2205 subroutine zsetmobius(z0,infz0, z1,infz1, z2,infz2, aa,bb,cc,dd)
2206 implicit none
2207 complex(dp), intent(IN ):: z0,z1,z2
2208 logical, intent(IN ):: infz0,infz1,infz2
2209 complex(dpc), intent(OUT):: aa,bb,cc,dd
2210 real(dp),parameter:: zero=0_dp,one=1_dp
2211 complex(dpc) :: z02,z10,z21
2212 z21=z2-z1
2213 z02=z0-z2
2214 z10=z1-z0
2215 if( (z0==z1.and.infz0.eqv.infz1).or.&
2216  (z1==z2.and.infz1.eqv.infz2).or.&
2217  (z2==z0.and.infz2.eqv.infz0)) &
2218  stop 'In setmobius; anchor points must be distinct'
2219 if(infz2 .or. (.not.infz0 .and. abs(z0)<abs(z2)))then
2220 ! z0 is finite and smaller than z2:
2221  if(infz1)then
2222  aa=one/sqrt(z02) ! <- z1 is infinite
2223  elseif(infz2)then
2224  aa=one/sqrt(z10) ! <- z2 is infinite
2225  else
2226  aa=sqrt(-z21/(z02*z10)) ! <- all zs are finite
2227  endif
2228  bb=-z0*aa
2229  if(infz1)then
2230  cc=aa ! <- z1 is infinite
2231  dd=-z2*aa !
2232  elseif(infz2)then
2233  cc=zero ! <- z2 is infinite
2234  dd=z10*aa !
2235  else
2236  cc=-(z10/z21)*aa ! <- all zs are finite
2237  dd= z2*(z10/z21)*aa !
2238  endif
2239 else
2240 ! z2 is finite and smaller than z0:
2241  if(infz0)then
2242  cc=one/sqrt(z21) ! <- z0 is inifinite
2243  elseif(infz1)then
2244  cc=one/sqrt(z02) ! <- z1 is infinite
2245  else
2246  cc=sqrt(-z10/(z02*z21)) ! <- all zs are finite
2247  endif
2248  dd=-z2*cc
2249  if(infz0)then
2250  aa=zero ! <- z0 is inifinite
2251  bb=-z21*cc !
2252  elseif(infz1)then
2253  aa=cc ! <- z1 is infinite
2254  bb=-z0*cc !
2255  else
2256  aa=(-z21/z10)*cc ! <- all zs are finite
2257  bb=z0*(z21/z10)*cc !
2258  endif
2259 endif
2260 end subroutine zsetmobius
2261 
2276 subroutine zmobius(aa,bb,cc,dd, z,infz, w,infw)! [mobius]
2277 implicit none
2278 complex(dpc),intent(IN ):: aa,bb,cc,dd,z
2279 logical, intent(IN ):: infz
2280 complex(dpc),intent(OUT):: w
2281 logical, intent(OUT):: infw
2282 real(dp),parameter:: zero=0_dp
2283 complex(dpc) :: top,bot
2284 w=0
2285 infw=.false.
2286 if(infz)then
2287  top=aa
2288  bot=cc
2289 else
2290  top=aa*z+bb
2291  bot=cc*z+dd
2292 endif
2293 
2294 if(abs(bot)==zero)then
2295  infw=.true.
2296 else
2297  w=top/bot
2298 endif
2299 end subroutine zmobius
2300 
2311 subroutine cmobius(aa,bb,cc,dd, vz,vw)! [mobius]
2312 implicit none
2313 complex(dpc), intent(IN ):: aa,bb,cc,dd
2314 real(dp),dimension(3),intent(IN ):: vz
2315 real(dp),dimension(3),intent(OUT):: vw
2316 complex(dpc):: z,w
2317 logical :: infz,infw
2318 call ctoz(vz, z,infz)
2319 call zmobius(aa,bb,cc,dd, z,infz, w,infw)
2320 call ztoc(w,infw, vw)
2321 end subroutine cmobius
2322 
2335 subroutine zmobiusi(aa,bb,cc,dd, zz,infz, zw,infw) ! [mobiusi]
2336 implicit none
2337 complex(dpc),intent(IN ):: aa,bb,cc,dd,zz
2338 logical, intent(IN ):: infz
2339 complex(dpc),intent(OUT):: zw
2340 logical, intent(OUT):: infw
2341 real(dp),parameter:: one=1_dp
2342 complex(dpc) :: aai,bbi,cci,ddi,d
2343 d=one/(aa*dd-bb*cc)
2344 aai=dd*d
2345 ddi=aa*d
2346 bbi=-bb*d
2347 cci=-cc*d
2348 call zmobius(aai,bbi,cci,ddi, zz,infz, zw,infw)
2349 end subroutine zmobiusi
2350 
2351 
2352 end module pmat4
real(sp) function, dimension(size(a), size(b)) outer_product_s(a, b)
Return the outer product matrix of two single precision real vectors.
Definition: pmat4.f90:272
real(sp) function huarea_s(sa, sb)
Spherical area of right-angle triangle whose orthogonal sides have orthographic projection dimensions...
Definition: pmat4.f90:629
subroutine dlltoxy_d(rlat, drlat, drlon, x2)
From a reference latitude, and increments of latitude and longitude, return the local cartesian 2-vec...
Definition: pmat4.f90:935
real(dp) function, dimension(size(a)) normalized_d(a)
Return the normalized version of a double precision real vector.
Definition: pmat4.f90:156
real(dp) function, dimension(3) axial33_d(b)
Return the 3-vector corresponding to the given antisymmetric &quot;axial vector&quot; matrix, assuming a right-handed correspondence.
Definition: pmat4.f90:461
real(dp) function hav_d(t)
Haversine function for double precision real (for geometry on the sphere).
Definition: pfun.f90:102
real(dp) function sarea_d(v1, v2, v3)
Compute the area of the spherical triangle, {v1,v2,v3}.
Definition: pmat4.f90:722
real(sp) function dtarea_s(rlat, drlata, drlona, drlatb, drlonb)
Compute the area of the spherical triangle with a vertex at latitude rlat, and two other vertices...
Definition: pmat4.f90:785
real(dp) function, dimension(size(b, 1)) diagnn_d(b)
Return the vector whose elements are the diagonal ones of a given matrix.
Definition: pmat4.f90:537
real(sp) function triple_product_s(a, b, c)
Return the triple product of three single precision real 3-vectors.
Definition: pmat4.f90:321
real(sp) function dqarea_s(rlat, drlata, drlona, drlatb, drlonb, drlatc, drlonc)
Compute the area of the spherical quadrilateral with a vertex at latitude rlat, and three other verti...
Definition: pmat4.f90:867
real(dp) function trace_d(b)
Return the trace of a given double precision real matrix.
Definition: pmat4.f90:578
real(sp) function trace_s(b)
Return the trace of a given single precision real matrix.
Definition: pmat4.f90:566
real(dp) function det_d(a)
Return the determinant of a double precision matrix.
Definition: pmat4.f90:368
integer(spi) function, dimension(size(a), size(a)) diagn_i(a)
Return the diagonal matrix whose elements are the given vector.
Definition: pmat4.f90:507
real(sp) function hav_s(t)
Haversine function for single precision real (for geometry on the sphere).
Definition: pfun.f90:89
subroutine expmatd(n, a, b, bd, detb, detbd)
Like expmat, but for the 1st derivatives also.
Definition: pmat4.f90:1766
real(dp) function dqarea_d(rlat, drlata, drlona, drlatb, drlonb, drlatc, drlonc)
Compute the area of the spherical quadrilateral with a vertex at latitude rlat, and three other verti...
Definition: pmat4.f90:896
real(sp) function, dimension(3) axial33_s(b)
Return the 3-vector corresponding to the given antisymmetric &quot;axial vector&quot; matrix, assuming a right-handed correspondence.
Definition: pmat4.f90:447
real(dp) function, dimension(size(u)) orthogonalized_d(u, a)
Return the part of vector a that is orthogonal to unit vector u.
Definition: pmat4.f90:186
real(sp) function absv_s(a)
Return the absolute magnitude of a single precision real vector.
Definition: pmat4.f90:118
real(sp) function, dimension(size(a)) normalized_s(a)
Return the normalized version of a single precision real vector.
Definition: pmat4.f90:142
subroutine gram_d(as, b, nrank, det)
Apply a form of Gram-Schmidt orthogonalization process to return as many normalized orthogonal basis ...
Definition: pmat4.f90:1077
subroutine normalize_s(v)
Normalize the given single precision real vector.
Definition: pmat4.f90:975
subroutine zmobiusi(aa, bb, cc, dd, zz, infz, zw, infw)
Perform the inverse of the mobius transformation with coefficients, {aa,bb,cc,dd}.
Definition: pmat4.f90:2335
real(dp) function, dimension(size(a), size(b)) outer_product_d(a, b)
Return the outer product matrix of two double precision real vectors.
Definition: pmat4.f90:288
integer(spi) function, dimension(size(a), size(b)) outer_product_i(a, b)
Return the outer product matrix of two integer vectors.
Definition: pmat4.f90:304
subroutine graml_d(as, b, nrank, detsign, ldet)
A version of gram_d where the determinant information is returned in logarithmic form (to avoid overf...
Definition: pmat4.f90:1149
real(sp) function, dimension(size(u)) orthogonalized_s(u, a)
Return the part of vector a that is orthogonal to unit vector u.
Definition: pmat4.f90:171
subroutine ztocd(z, infz, v, vd)
The convention adopted for the complex derivative is that, for a complex infinitesimal map displaceme...
Definition: pmat4.f90:2082
subroutine plaingram_d(b, nrank)
A &quot;plain&quot; (unpivoted) version of Gram-Schmidt, for square matrices only.
Definition: pmat4.f90:1262
subroutine normalize_d(v)
Normalize the given double precision real vector.
Definition: pmat4.f90:987
integer(dpi) function det_id(a)
Return the determinant of a double precision integer matrix.
Definition: pmat4.f90:403
Module for handy vector and matrix operations in Euclidean geometry.
Definition: pmat4.f90:45
real(sp) function det_s(a)
Return the determinant of a single precision matrix.
Definition: pmat4.f90:347
real(sp) function, dimension(size(b, 1)) diagnn_s(b)
Return the vector whose elements are the diagonal ones of a given matrix.
Definition: pmat4.f90:522
Standard integer, real, and complex single and double precision kinds.
Definition: pkind.f90:7
real(sp) function, dimension(3, 3) axial3_s(a)
Return the axial &quot;vector&quot;, as an antisymmetric matrix, corresponding to the given 3-vector assuming a...
Definition: pmat4.f90:419
real(dp) function triple_product_d(a, b, c)
Return the triple product of three double precision real 3-vectors.
Definition: pmat4.f90:335
subroutine cmobius(aa, bb, cc, dd, vz, vw)
Perform a complex Mobius transformation from cartesian vz to cartesian vw where the transformation co...
Definition: pmat4.f90:2311
subroutine expmatdd(n, a, b, bd, bdd, detb, detbd, detbdd)
Like expmat, but for the 1st and 2nd derivatives also.
Definition: pmat4.f90:1833
real(sp) function, dimension(size(a), size(a)) diagn_s(a)
Return the diagonal matrix whose elements are the given vector.
Definition: pmat4.f90:475
integer(spi) function trace_i(b)
Return the trace of a given integer matrix.
Definition: pmat4.f90:590
subroutine zsetmobius(z0, infz0, z1, infz1, z2, infz2, aa, bb, cc, dd)
Find the Mobius transformation complex coefficients, aa,bb,cc,dd, with aa*dd-bb*cc=1, that takes polar stereographic point, z0 to the north pole, z1 to (lat=0,lon=0), z2 to the south pole (=complex infinity).
Definition: pmat4.f90:2205
Some of the commonly used constants (pi etc) mainly for double-precision subroutines.
Definition: pietc.f90:14
real(dp) function absv_d(a)
Return the absolute magnitude of a double precision real vector.
Definition: pmat4.f90:130
subroutine rowgram(m, n, a, ipiv, tt, b, rank)
Without changing (tall) rectangular input matrix a, perform pivoted gram- Schmidt operations to ortho...
Definition: pmat4.f90:1314
real(dp) function, dimension(3, 3) axial3_d(a)
Return the axial &quot;vector&quot;, as an antisymmetric matrix, corresponding to the given 3-vector assuming a...
Definition: pmat4.f90:433
real(dp) function, dimension(4) triple_cross_product_d(u, v, w)
Return the triple_cross_product for 4-vectors.
Definition: pmat4.f90:252
Definition: pmat.f90:20
real(sp) function, dimension(3) cross_product_s(a, b)
Return the cross product of two single precision real 3-vectors.
Definition: pmat4.f90:201
subroutine plaingram_s(b, nrank)
A &quot;plain&quot; (unpivoted) version of Gram-Schmidt, for square matrices only.
Definition: pmat4.f90:1226
integer(spi) function det_i(a)
Return the determinant of a single precision integer matrix.
Definition: pmat4.f90:389
real(dp) function, dimension(size(a), size(a)) diagn_d(a)
Return the diagonal matrix whose elements are the given vector.
Definition: pmat4.f90:491
subroutine dlltoxy_s(rlat, drlat, drlon, x2)
From a reference latitude, and increments of latitude and longitude, return the local cartesian 2-vec...
Definition: pmat4.f90:915
integer(spi) function, dimension(3, 3) identity3_i()
Return the 3-dimensional integer identity matrix.
Definition: pmat4.f90:614
integer(spi) function, dimension(n, n) identity_i(n)
Return the integer identity matrix for a given dimensionality.
Definition: pmat4.f90:602
integer(spi) function, dimension(size(b, 1)) diagnn_i(b)
Return the vector whose elements are the diagonal ones of a given matrix.
Definition: pmat4.f90:552
real(sp) function, dimension(4) triple_cross_product_s(u, v, w)
Deliver the triple-cross-product, x, of the three 4-vectors, u, v, w, with the sign convention that o...
Definition: pmat4.f90:231
real(dp) function huarea_d(sa, sb)
Spherical area of right-angle triangle whose orthogonal sides have orthographic projection dimensions...
Definition: pmat4.f90:647
subroutine gram_s(as, b, nrank, det)
Apply a form of Gram-Schmidt orthogonalization process to return as many normalized orthogonal basis ...
Definition: pmat4.f90:1006
real(sp) function sarea_s(v1, v2, v3)
Compute the area of the spherical triangle, {v1,v2,v3}, measured in the right-handed sense...
Definition: pmat4.f90:667
subroutine zmobius(aa, bb, cc, dd, z, infz, w, infw)
Perform a complex Mobius transformation from (z,infz) to (w,infw) where the transformation coefficien...
Definition: pmat4.f90:2276
real(dp) function dtarea_d(rlat, drlata, drlona, drlatb, drlonb)
Compute the area of the spherical triangle with a vertex at latitude rlat, and two other vertices...
Definition: pmat4.f90:824
subroutine cyclic(u1, u2, u3, d1, d2, d3)
Cyclically permute real vectors, u1, u2, u3, and scalars, d1, d2, d3.
Definition: pmat4.f90:704
real(dp) function, dimension(3) cross_product_d(a, b)
Return the cross product of two double precision real 3-vectors.
Definition: pmat4.f90:214