!SA2-hVecMath16Mod.f95
!2025.05.24.2030cdt- Dodecahedron Analysis & other real(16) subroutines.

!            Author- Jeffrey M. Setterholm, Lakeville,MN 55044 USA
!         IP Status- Free source code (e.g.: post copyright)
!
!          Computer- "T3"/Dell Precision T3500/Intel i5 E5520/Win10Pro-21H2
!                     ^name ^Mfgr.Id           ^chipset       ^OS 
!                    /Absoft Pro Fortran 21.0.2/GeForce GTX 1050/f90gl~Glut3.7
!                     ^compiler ~Fortran 95     ^graphics card   ^graphics

!    f90gl bindings- public domain; see "https://math.nist.gov/f90gl/"

!Disclaimer:
! ***************************************************************
! **********   Individual cognition is always flawed,  **********
! **********         including yours and mine.         **********
! **********                 - So: -                   **********
! **********       Use this code at your own risk.     **********
! ***************************************************************

!Table of Contents: ...use to search...
!-----
!Module hVecMath16Mod
! Module DodecIO
! contains
!  Function Sqrj(Valin) Result(ValSqrout)                       2025.04.21.1545
!End Module DodecIO
!Subroutine DodecModel16(iP)                                    2013.10.23.0840
!Subroutine SinCos16(AngleD,SinAngle,CosAngle)                  2013.10.22.0940
!Subroutine ASinD16(Angle,SineAngle,iP)                         2013.10.20.0810
!Subroutine ACosD16(Angle,CosineAngle,iP)                       2013.10.20.0810
!Subroutine hPointPolar16(Ph,aRpyh,PDmag)                       2013.10.17.2210
!Subroutine hVnorm16(Vinh,Vouth,Vmag)                           2013.10.17.2210
!Subroutine hVecMirror2D16(Vinh,Vmirrorh,Vouth)                 2013.10.09.0800
!Subroutine hVecMirror3D16(Vinh,Vmirror1h,Vmirror2h,Vouth)      2013.10.09.0800
!Subroutine hCross16(Voh,Vi1h,Vi2h)                             2013.10.09.1005
!Subroutine MantExp16(ValIn16,rMant16,iExp,iP)                  2025.04.27.1110
!Subroutine Sqrt16a(ValIn,ValSqrOut,iMZP,iP)      Rev.A         2025.04.27.1110
!Function Factrl(N)                                             2025.04.21.1510
!------------------------------------------------------------------------------
!---------------------------------------------------------------------------7-9
!------------------------------------------------------------------------------

Module DodecIO    !This is the user interface to the Dodecahedron model.
!2013.10.23.0840cdt JMS- Dodecahedron Input/Output Variables
!                      - DodecIO: Now using "type(DodecRec),save::Dd"
!use DodecIO
!use DodecIO, only:vNormal,Vertex,iFace,iEdge,DdWicki,iWikiToSett

 type::  DodecRec  ;sequence                                        !2013.10.16
  integer(4)::InitDodec
  !Wikipedia's Parameters - Theoretical values:
  real(16)  ::EdgeL          !Edge length   =(sqrt(5)-1)/sqrt(3) exactly
                            !~= 0.713644179546179863883939686092_16
  !For an edge of length 1   Radii:
  real(16)  ::ru             !outer           = (1 +  sqrt(5))  *sqrt(3) /4
  real(16)  ::rm             !mid-edge        = (     sqrt(5)      +3   )/4
  real(16)  ::ri             !inner           = (sqrt(sqrt(5)*11/10+5/2))/2
  real(16)  ::P1             !of face vertices=  sqrt(ru*ru-ri*ri)
  !Wikipedia's vertices are on a sphere with radius=sqrt(3):    
  real(16)  ::Phi            ! =(1+sqrt(5))/2 exactly = the "Golden Ratio"
                            ! = 1.61803398874989484820458683436564_16
  real(16)  ::DPhi           ! ="one divided by Phi"
                            !~= 0.61803398874989484820458683436564_16
  real(16)  ::DdWicki(4,20)
  !...the resulting coordinates must be divided by sqrt(3)
  !...and pitched down:
  real(16)  ::WickiPitch    !=90-(acosd(rm/ru)+acosd(ri/ru)
                           !=-31.71747441146100532421390313982794_16 degrees
                           !                            ...982700 ...playback
  real(16)  ::EdgeLengthWiki(30) !Edge- lengths- ~= EdgeL ...Yes!
  real(16)  ::EdgeLErrorRmsWiki  !           - RMS error
  !...and renumbered as follows:
  integer(4)::iWikiToSett(20) 
                  !=(/19,16,4,2,15,7,14,18,13,5,11,3,17,1,9,20,8,12,10,6/)
  !... to match my model.    

  !So, my model parameters:
                            !... are theoretically quantified by:
  real(16)  :: sinV1Ang      != (ri-2*P1)/(ru*sqrt(5))
  real(16)  ::cosV16Ang      !=       P1 / ru
  real(16)  :: v1Ang         !Vertex -  #1- Pitch Angle w.r.t. the X-Y plane.
    !V1Ang   =-10.81231696357170629128494477447353_16 !...Theoretical
    !V1Ang   =-10.81231696357170629128494477447422_16 !...best fit to date.
  real(16)  ::v16Ang         !       - #16-
    !V16Ang  =-52.62263185935030435714286150510032_16 !...Theoretical
    !V16Ang  =-52.62263185935030435714286150501010_16 !...best fit to date.
  !Indices:
  integer(4)::iFace(5,12)    !Face   - vertex  # list (CCW looking out)
  integer(4)::iEdge(2,30)    !Edge   - vertex  # list
  integer(4)::iVertex(3,20)  !vertex - vNormal # list (CCW looking out)
  integer(4)::iEdgevN(2,30)  !edge   - vNormal # list
  !Computed Vertices & Normals:
  real(16)  ::Vertex(4,20)   !   Vertex   vector- (X,Y,Z,1.)
  real(16)  ::vNormal(4,12)  !Face normal vector- (X,Y,Z,1.)
  real(16)  ::EdgeLength(30) !Edge- lengths- ~= EdgeL ...Yes!
  real(16)  ::EdgeVnDist(30) !Edge- Distance between vNormal vectors
  real(16)  ::EdgeLErrorRms  !             - RMS error
  !Nominal number of significant digits for real(16) variables herein: 
  integer(4)::iQuadDigits    !=31 for AbsoftProFortran vsn 9.0 (used here).
  !Graphics drawing flags: 
  integer(4)::  iFaceSee(0:12)=1 !if(iFaceSee(0)=0: set desired face flags=1
                            !=(/0,1,0,0,0,0,1,0,0,0,0,1,1/) for #1,#6,#11,&#12
  integer(4)::iVertexSee(0:20)=0 !computed- based on iFaceSee(1:12).
  integer(4)::  iEdgeSee(0:30)=0 !        -
                      ! ^      (0)=-1: disable labels, (0)=-2: disable lines
 end type DodecRec ;type(DodecRec),save::Dd !<-for "Dodecahedron"
!---------------------------------------------------------------------------7 9
 Contains
!---------------------------------------------------------------------------7 9

Function Sqrj(Valin) Result(ValSqrout) !<-real(16)
!2025.04.21.1545cdt JMS- more accurate than intrinsic function sqrt(Valin).
!                      - ~ 30-place accuracy.
implicit none                                                        !arguments
  real(16)::ValIn  !Input number
  real(16)::ValSqrout != square root of Valin
 !--                                                                 !internals
  integer(4)::iP
  integer(4) iMZP
  !-----
                               iP = 0
   call Sqrt16a(ValIn,ValSqrOut,iMZP,iP)
End Function Sqrj
!---------------------------------------------------------------------------7 9

End Module DodecIO
!------------------------------------------------------------------------------
!---------------------------------------------------------------------------7 9
!------------------------------------------------------------------------------
!Module hVecMath16Mod
!use hVecMath16Mod                                                  !2011.04.25
!---
!Implicit none
!---
!public ::        & !*****************************************
!    DodecModel16   & !(iP)                                    *2013.10.18.1910
! ,  hPointPolar16  & !(Ph,aRpyh,PDmag)                        *2013.10.17.2210
! ,  hVnorm16       & !(Vinh,Vouth,Vmag)                       *2013.10.17.2210
! ,  hVecMirror2D16 & !(Vinh,Vmirrorh,Vouth)                    2013.10.09.0800
! ,  hVecMirror3D16 & !(Vinh,Vmirror1h,Vmirror2h,Vouth)         2013.10.09.0800
! ,  SinCos16       & !(AngleD,SinAngle,CosAngle)               2013.10.14.1150
! ,  ASinD16        & !(Angle,SineAngle,iP)                    *2013.10.18.1315
! ,  ACosD16        & !(Angle,CosineAngle,iP)                  *2013.10.18.1315
! ,  hCross16       & !(Voh,Vi1h,Vi2h)                          2013.10.09.1005
! ,  Sqrt16           !(ValIn,ValSqrOut,iP)                    *2013.10.18.1000

!---------------------------------------------------------------------------7-9

Subroutine DodecModel16(iP)
!2013.10.22.1950cdt JMS- Graphics drawing flags.
!2013.10.20.1405cdt JMS- Using best-fit value of V1Ang & V16Ang
!2013.10.16.1815cdt JMS- Dodecahedron - the 12-face Platonic Solid
!                         which fits inside a sphere of radius=1.0
!------                - Reference: Wikipedia's article on "Dodecahedron".
 use DodecIO, only:Dd,Sqrj !,Pi16
 implicit none
 !-----                                                               Arguments
  integer(4)::iP
 !-----                                                               Internals
  !-- Internal use variables
  integer(4)::i,iDel,iDel2,j,k,L
  real(16)  ::v0=0._16  !0.
  real(16)  ::v1=1._16  !1.
  real(16)  ::Angle     !Yaw angle-
  real(16)  ::Ver1(4)   !Vertex vector- # 1- used as a template
  real(16)  ::Ver16(4)  !             - #16-
  real(16)  ::sY,cY     !Yaw  - cosine, sine
  real(16)  ::sP,cP     !Pitch- cosine, sine
  real(16)  ::vTemph(4) !Temporary vector
  real(16)  ::DdWickiTemp(4,20)
  real(16)  ::Angle1,Angle2
  integer(4)::iReverse(20)
  integer(4)::iMZP

  real(16)::Vsqrt                                                   !2025.04.20
  real(16)::Pi                                                      !2025.04.20

   !----- Dodecahedron initilization - Jeff Setterholm's algorithm:
   if(Dd%InitDodec.eq.0) then !-------------------------------
     if(iP>5) write(iP, &
    "('Dodecahedron model initialization: SA2-hVecMath16Mod.f95 @L181+')")
     !---
     Dd%iQuadDigits=precision(Dd%EdgeL)
     !---
     Dd%EdgeL  =(sqrj(5._16)-1._16)/sqrj(3._16)       !Theoretical value.
      !=   0.71364417954617986388393968609217_16

     if(0>1) then !---------- Bypasses when (0>1); enters using (1>0)
       if(iP>5) then  !Testing Sqrt(), call sqrt16(), & sqrj()
       write(iP,"(/,'Testing Absoft`s sqrt():')")
       Vsqrt   = sqrt(5._16)
       write(iP,"(e45.35,' :5.-sqrt(  5.)*sqrt(  5.)'  )") 5._16-Vsqrt*Vsqrt
      !Testing Absoft`s sqrt():
      !   -0.40245584642661924590356647968292230E-15 :5.-sqrt(  5.)*sqrt(  5.)
       write(iP,"(  'Testing my sqrt16() subroutine:')")
       call sqrt16a(5._16,Vsqrt,iMZP,0)
       write(iP,"(e45.35,' :5.-sqrt16(5.)*sqrt16(5.)'  )") 5._16-Vsqrt*Vsqrt
      !Testing my sqrt16() subroutine:
      !    0.00000000000000000000000000000000000E+00 :5.-sqrt16(5.)*sqrt16(5.)
       write(iP,"(  'Testing my sqrj() function:')")
       Vsqrt = sqrj(5._16)
       write(iP,"(e45.35,' :5.-sqrj(  5.)*sqrj(  5.)',/)") 5._16-Vsqrt*Vsqrt
      !Testing my sqrj() function:
      !    0.00000000000000000000000000000000000E+00 :5.-sqrj(  5.)*sqrj(  5.)
       write(iP,"('Dd%EdgeL = '/f40.35)")  Dd%EdgeL
       endif!(iP>5)
     endif!(0>1) -------- end bypass ------------------------------------------


     Dd%ru= (1._16+sqrj(5._16))*sqrj(3._16)/4._16
      !=  1.40125853844407354467667793532207_16
     Dd%rm= (sqrj(5._16)+3._16)/4._16
      !=  1.30901699437494742410229341718282_16
     Dd%ri= (sqrj(sqrj(5._16)*11._16/10._16+5._16/2._16))/2._16
      !=  1.11351636441160673519437503948696_16 
     Dd%P1=  sqrj(Dd%ru*Dd%ru-Dd%ri*Dd%ri)
      !=  
    !-- wikipedia's vertex coordinates:
     Dd%Phi     = (1._16+sqrj(5._16))/2._16
      !=  1.61803398874989484820458683436564_16
     Dd%dPhi=  1._16/Dd%Phi
      !=  0.61803398874989484820458683436564_16
     call ACosD16(Angle1,(Dd%rm/Dd%ru),0)
     call ACosD16(Angle2,(Dd%ri/Dd%ru),0)
     Dd%WickiPitch=90._16-(Angle1+angle2)
      !=-31.71747441146100532421390313982794_16

    !Dd%sinV1Ang=(   ri-2    *   P1)/(   ru*sqrt(5    ))
     Dd%sinV1Ang=(Dd%ri-2._16*Dd%P1)/(Dd%ru*sqrj(5._16))
               != -0.18759247408507989986013934690761_16
     call ASinD16(Dd%V1Ang,Dd%sinV1Ang,0)
    !Dd%V1Ang   =-10.81231696357170629128494477447353_16 !=Pitch ang.v#1 theor.
     Dd%V1Ang   =-10.81231696357170629128494477447422_16 !...best fit to date.
    !            -10.81231696357170629128494477447412    !...as reported back
    !                                          ^ 32nd place 
     Dd%cosV16Ang=Dd%P1/Dd%ru
     call ACosD16(Dd%V16Ang,Dd%cosV16Ang,0) ;
                  Dd%V16Ang=-Dd%v16Ang    
    !Dd%V16Ang  =-52.62263185935030435714286150510032_16 !=Pitch ang.v#16 theor
     Dd%V16Ang  =-52.62263185935030435714286150501010_16 !...best fit to date.
    !            -52.62263185935030435714286150500763_16 !...as reported back
    !                                          ^ 32nd place 

    !-- Index faces vertices, edge vertices, vertex faces, & edge faces:
     iDel =5 ;iDel2=0
     do i=1,5 !Index faces(10 of 12 sets) {-CCW looking out} & all edges(30):
         Dd%iFace(1:5,i    )=(/i  ,i+5,i+10,i+ 9+iDel ,i+4+iDel /) ![ 1, 5]
         Dd%iFace(1:5,i  +5)=(/i+5,i  ,i+15,i+16-iDel2,i+1-iDel2/) ![ 6,10]
         Dd%iEdge(1:2,i*6-5)=  Dd%iFace(1:2,i)                   !1, 7,13,19,25
         Dd%iEdge(1:2,i*6-4)=  Dd%iFace(2:3,i)                   !2, 8,14,20,26
         Dd%iEdge(1:2,i*6-3)=  Dd%iFace(3:4,i)                   !3, 9,15,21,27
         Dd%iEdge(1:2,i*6-2)=(/Dd%iFace(1,i+5),Dd%iFace(5,i+5)/) !4,10,16,22,28
         Dd%iEdge(1:2,i*6-1)=(/Dd%iFace(5,i+5),Dd%iFace(4,i+5)/) !5,11,17,23,29
         Dd%iEdge(1:2,i*6  )=(/Dd%iFace(4,i+5),Dd%iFace(3,i+5)/) !6,12,18,24,30
       Dd%iVertex(1:3,i    )=(/i  ,i+4+iDel ,i+5/)               ![ 1, 5]
       Dd%iVertex(1:3,i  +5)=(/i+5,i+1-iDel2,i  /)               ![ 6,10]
       Dd%iVertex(1:3,i +10)=(/i  ,i+1-iDel2, 11/)               ![11,15]
       Dd%iVertex(1:3,i +15)=(/i+5,i+4+iDel , 12/)               ![16,20]

       Dd%iEdgevN(1:2,i*6-5)=(/i        ,i+ 5      /)            !1, 7,13,19,25
       Dd%iEdgevN(1:2,i*6-4)=(/i        ,i+ 1-iDel2/)            !2, 8,14,20,26
       Dd%iEdgevN(1:2,i*6-3)=(/i        ,  11      /)            !3, 9,15,21,27
       Dd%iEdgevN(1:2,i*6-2)=(/i+1-iDel2, i+5      /)            !4,10,16,22,28
       Dd%iEdgevN(1:2,i*6-1)=(/i+5      , i+6-iDel2/)            !5,11,17,23,29
       Dd%iEdgevN(1:2,i*6  )=(/i+5      ,  12      /)            !6,12,18,24,30
       if(i.eq.1) iDel =0
       if(i.eq.4) iDel2=5
     enddo!i
         Dd%iFace(  1:5,11)=(/13,14,15,11,12/)                 !Face #11-Top
         Dd%iFace(  1:5,12)=(/19,18,17,16,20/)                 !Face #12-Bottom

     !-- Wikipedia model: ---------------------------------
     Dd%iWikiToSett=(/19,16,4,2,15,7,14,18,13,5,11,3,17,1,9,20,8,12,10,6/)
                                                 L=0
     do i=-1,1,2 ;do j=-1,1,2 ;do k=-1,1,2      ;L=L+1
     Dd%DdWicki(1:4,L)=(/i*1._16  ,j*1._16  ,k*1._16  ,1._16/)
     enddo       ;enddo       ;enddo
     do j=-1,1,2 ;do k=-1,1,2                   ;L=L+1
     Dd%DdWicki(1:4,L)=(/  0._16  ,j*Dd%DPhi,k*Dd%Phi ,1._16/)
     enddo       ;enddo
     do i=-1,1,2 ;do j=-1,1,2                   ;L=L+1
     Dd%DdWicki(1:4,L)=(/i*Dd%DPhi,j*Dd%Phi ,0._16    ,1._16/)
     enddo       ;enddo
     do i=-1,1,2              ;do k=-1,1,2      ;L=L+1
     Dd%DdWicki(1:4,L)=(/i*Dd%Phi ,  0._16  ,k*Dd%DPhi,1._16/)
     enddo                    ;enddo
     !Normalize the vertex vectors:
     do L=1,20 ;Dd%DdWicki(1:3,L)=Dd%DdWicki(1:3,L)/sqrj(3._16) ;enddo!L

   !-- Compute the Edge Lengths for Wikipedia's model:

     !Rotate the vertex vectors to the "Setterholm Model" orientation:
     call SinCos16(Dd%WickiPitch,sp,Cp)
     do j=1,20 ;L=Dd%iWikiToSett(j) ; iReverse(L)=j
       vTemph(1)= Dd%DdWicki(1,L)*cP+Dd%DdWicki(3,L)*sP 
       vTemph(2)= Dd%DdWicki(2,L)
       vTemph(3)=-Dd%DdWicki(1,L)*sP+Dd%DdWicki(3,L)*cP
       DdWickiTemp(1:3,j)=vTemph(1:3)
     enddo!L
     !Remap the Wickimodel vector indices.
     Dd%DdWicki=DdWickiTemp

    !if(iP>5) write(iP,"('Edge length errors @L303:')")

     Dd%EdgeLErrorRmsWiki=0._16
     do i=1,30   
       vTemph= Dd%DdWicki(1:4,Dd%iEdge(2,i))  &
              -Dd%DdWicki(1:4,Dd%iEdge(1,i))
       Dd%EdgeLengthWiki(i)= ( vTemph(1)*vTemph(1) &
                              +vTemph(2)*vTemph(2) &
                              +vTemph(3)*vTemph(3)   )
       Dd%EdgeLengthWiki(i)=sqrj(Dd%EdgeLengthWiki(i))-Dd%EdgeL
       Dd%EdgeLErrorRmsWiki=Dd%EdgeLErrorRmsWiki+ Dd%EdgeLengthWiki(i)  &
                                                 *Dd%EdgeLengthWiki(i)

      !write(iP,"(i3,f40.35)") i,Dd%EdgeLErrorRmsWiki !/i
     enddo!i
     Dd%EdgeLErrorRmsWiki=sqrj(Dd%EdgeLErrorRmsWiki/30._16)
    !write(iP,"(/3x,f40.35,' Dd%EdgeLErrorRmsWiki')")  Dd%EdgeLErrorRmsWiki


     Dd%InitDodec=1 ;write(6,"(a1)") char(7)                              !Beep
   endif !InitDodec=0 ---------------------------------------------------------

   !-- Update the Graphics drawing flags:
  !Dd%iFaceSee(0:12)=(/0,1,0,0,0,0,1,0,0,0,0,1,1/)
   Dd%iVertexSee(1:20)=0 ;Dd%iEdgeSee(1:30)=0
   if(Dd%iFaceSee(0).eq.1) Dd%iFaceSee=1                                 !Reset
   do i=1,12 ;if(Dd%iFaceSee(i).eq.0) cycle
     do j=1,5 ;k=Dd%iFace(j,i) ;Dd%iVertexSee(k)=1 ;enddo!j
   enddo!i
   do i=1,30 ;if(Dd%iVertexSee(Dd%iEdge(1,i)).eq.0) cycle
              if(Dd%iVertexSee(Dd%iEdge(2,i)).eq.0) cycle
     Dd%iEdgeSee(i)=1
   enddo!i

   !----- From here on: V1Ang & V16Ang quantify the model...
   !-- Setup the Vertex templates:
   call SinCos16(Dd%v1Ang ,sP,cP) ;Ver1 =(/cP,v0,-sP,v1/) != vertex#1
   call SinCos16(Dd%v16Ang,sP,cP) ;Ver16=(/cP,v0,-sP,v1/) != vertex#16
   do i=1,5  !Spin the vertex templates to define all the vertices:
     Angle=(i-1)*72._16       ;call SinCos16(Angle,sY,cY)
     Dd%Vertex( 1:4,i  )=(/ Ver1(1)*cY, Ver1(1)*sY,  Ver1(3),v1/)! 1- 5
     Dd%Vertex(1:4,i+15)=(/Ver16(1)*cY,Ver16(1)*sY, Ver16(3),v1/)!16-20
     Angle=(i-1)*72._16+36._16 ;call SinCos16(Angle,sY,cY)
     Dd%Vertex( 1:4,i+5)=(/ Ver1(1)*cY, Ver1(1)*sY, -Ver1(3),v1/)! 6-10
     Dd%Vertex(1:4,i+10)=(/Ver16(1)*cY,Ver16(1)*sY,-Ver16(3),v1/)!11-15
   enddo!i

   !-- Compute each Face Normal vector by averaging its five vertices:
   do i=1,12 ;Dd%vNormal(1:4,i)=0._16
     do j=1,5
       Dd%vNormal(1:4,i)=Dd%vNormal(1:4,i)+Dd%Vertex(1:4,Dd%iFace(j,i))
     enddo!j
     Dd%vNormal(1:4,i)=Dd%vNormal(1:4,i)/5._16
   enddo!i

   !-- Compute the resulting Edge Lengths:
   Dd%EdgeLErrorRms=0._16
   do i=1,30
     vTemph=Dd%Vertex(1:4,Dd%iEdge(2,i))-Dd%Vertex(1:4,Dd%iEdge(1,i))
     Dd%EdgeLength(i)= ( vTemph(1)*vTemph(1) &
                        +vTemph(2)*vTemph(2) &
                        +vTemph(3)*vTemph(3)   )
     call Sqrt16a(Dd%EdgeLength(i),Dd%EdgeLength(i),iMZP,0)
     Dd%EdgeLErrorRms=Dd%EdgeLErrorRms+ (Dd%EdgeLength(i)-Dd%EdgeL)  &
                                       *(Dd%EdgeLength(i)-Dd%EdgeL)
   enddo!i
   Dd%EdgeLErrorRms=sqrj(Dd%EdgeLErrorRms/30._16)

   !----- Report the results:
   if(iP.gt.5) then
     write(iP,"(/79('/'))")
     write(iP,"( 'Begin DodecModel16 printout @L374+' )")
     write(iP,"(/'  real(16) precision =',i3,' significant digits')")  &
                                      Dd%iQuadDigits
     write(iP,"( '     & tentatively: printouts truncate at 34 places max.')")
     write(iP,"(/'ru            =',f40.32)")  Dd%ru
     write(iP,"( 'rm            =',f40.32)")  Dd%rm
     write(iP,"( 'ri            =',f40.32)")  Dd%ri
     write(iP,"( 'P1            =',f40.32)")  Dd%P1
     write(iP,"(/' sinV1Ang     =',f40.32,' -theoretical'    )")    Dd%sinV1Ang
     write(iP,"( ' V#1 pitch Ang=',f40.32,' -best-fit value used' )")  Dd%V1Ang
     write(iP,"( 'cosV16Ang     =',f40.32,' -theoretical'    )")   Dd%cosV16Ang
     write(iP,"( 'V#16 Pitch Ang=',f40.32,' -best-fit value used'/)") Dd%V16Ang
     write(iP,"( 'Vertices:')")
       write(iP,"(i2,' =(/',2(f35.31,','),f35.31,'/)')")   1, ver1( 1:3)
       write(iP,"(i2,' =(/',2(f35.31,','),f35.31,'/)'/)") 16, ver16(1:3)
     do i=1,20
       write(iP,"(i2,' =(/',2(f35.31,','),f35.31,'/)')") i,Dd%Vertex( 1:3,i)
     enddo!i
     write(iP,"(/'Face normals:')")
     do i=1,12
       write(iP,"(i2,' =(/',2(f35.31,','),f35.31,'/)')") i,Dd%vNormal( 1:3,i)
     enddo!i
     !-- Indices & Lengths:
     write(iP,"(/'       CCW looking out    CCW looking out' )")
     write(iP,"( 'Face:    VertexIndex   Vertex:  FaceIndex ',\)")
     write(iP,"( ' Edge:   vNIndex VertexIndex &   ComputedEdgeLength')")
     do i=1,30
       if(i.le.12) then
         write(iP,"(' #',i2,'=( ',4(i2,','),i2,'/)',\)") i,Dd%iFace(1:5,i)
        else
         write(iP,"('                       ',\)")
       endif !i<=12
       if(i.le.20) then
         write(iP,"('   #',i2,'=( ',2(i2,','),i2,'/)',\)") i,Dd%iVertex(1:3,i)
        else
         write(iP,"('                   ',\)")
       endif !i<=20
       write(iP,"('   #',i2,' =( ',1(i2,','),i2,'/)',\)") i,Dd%iEdgeVn(1:2,i)
       write(iP,"(          ' =( ',1(i2,','),i2,'/)',\)")   Dd%iEdge(  1:2,i)
       write(iP,"(' :',f35.32)") Dd%EdgeLength(i)
     enddo!i
       write(iP,"(/48x,'Theoretical Edge Length=',f37.34)") Dd%EdgeL
       write(iP,"( 84x,'1         2         3')")
       write(iP,"( 75x,'12345678901234567890123456789012')")
       write(iP,"( 79x,'RMS Error=',e13.6,' : ^')") Dd%EdgeLErrorRms

     Dd%EdgeL  =(sqrj(5._16)-1._16)/sqrj(3._16)       !Theoretical value.
      !=   0.71364417954617986388393968609217_16

     write(iP,"( 'Wikipedia Model:')")
     write(iP,"( 'Phi           =',f40.37)")  Dd%Phi
     !=      1.61803398874989484820458683436564_16
     write(iP,"( 'dPhi          =',f40.37)")  Dd%dPhi
     !=      0.61803398874989484820458683436564_16
     write(iP,"( '...vertex algebraic summary:')")
     write(iP,"( 'W#:  S#:     Numeric values:')")
                                                 L=0
     do i=-1,1,2 ;do j=-1,1,2 ;do k=-1,1,2      ;L=L+1
       write(iP,"(i3,' =',i3,' =(/',i2,'     ,',i2,'     ,',i2,'     /)')")  &
       L,iReverse(L)        ,i          ,j          ,k
     enddo       ;enddo       ;enddo
     do j=-1,1,2 ;do k=-1,1,2                   ;L=L+1
       write(iP,"(i3,' =',i3,' =(/',i2,'     ,',i2,'*dPhi,',i2,'*phi /)')")  &
       L,iReverse(L)        ,0          ,j          ,k
     enddo       ;enddo
     do i=-1,1,2 ;do j=-1,1,2                   ;L=L+1
       write(iP,"(i3,' =',i3,' =(/',i2,'*dPhi,',i2,'*Phi ,',i2,'     /)')")  &
       L,iReverse(L)        ,i          ,j          ,0
     enddo       ;enddo
     do i=-1,1,2              ;do k=-1,1,2      ;L=L+1
       write(iP,"(i3,' =',i3,' =(/',i2,'*Phi ,',i2,'     ,',i2,'*dPhi/)')")  &
       L,iReverse(L)        ,i          ,0          ,k
     enddo                    ;enddo
      !write(iP,"( 37x,'RMS Edge Length Errors before rotation=',d13.6)") &
      !                 EdgeLErrorRmsWiki
     write(iP,"(/'...scale and rotate the model:')")
    !write(iP,"( 'WickiPitch    =',f40.37)")  WickiPitch
     !=    -31.71747441146100532421390313982794_16
     write(iP,"( 'Divide by sqrj(3)=',f40.32)")  sqrj(3._16)
     !=    -31.71747441146100532421390313982794_16
     write(iP,"( 'and pitch down by ',f40.32,' degrees.')")  Dd%WickiPitch
     !=    -31.71747441146100532421390313982794_16

     write(iP,"(/'... resulting Wicki-Vertices in Setterholm Model order:')")
     write(iP,"( 'W#:')")
     do i=1,20 ;L=Dd%iWikiToSett(i) !=iReverse(i)
       write(iP,"(i2,' =(/',2(f35.31,','),f35.31,'/)')") L,Dd%DdWicki( 1:3,i)
     enddo!i
     write(iP,"( '       Edge#   W1 W2       Length',29x,'Length error')")
     do i=1,30
       write(iP,"(9x,'#',i2,'=(/',i2,',',i2,'/)=',f38.34,d13.4)") &
       i,Dd%iWikiToSett(Dd%iEdge(1:2,i)),Dd%EdgeLengthWiki(i)+Dd%EdgeL &
                                        ,Dd%EdgeLengthWiki(i)
     enddo!i
       write(iP,"(     'Theoretical Edge Length=',f37.34)") Dd%EdgeL
       write(iP,"( 36x,'1         2         3')")
       write(iP,"( 27x,'1234567890123456789012345678901234')")
       write(iP,"( 33x,'RMS Error=',e13.6,' : ^')") Dd%EdgeLErrorRmsWiki

     write(iP,"(/'Graphics drawing flags:' )")
     write(iP,"( i2,1x,12i1,18x,' :  iFaceSee')")   Dd%iFaceSee
     write(iP,"( i2,1x,20i1,10x,' :iVertexSee')") Dd%iVertexSee
     write(iP,"( i2,1x,30i1    ,' :  iEdgeSee')")   Dd%iEdgeSee
     write(iP,"( ' ^          1         2         3')")
     write(iP,"( ' 0|123456789012345678901234567890')")
     write(iP,"(/'End of DodecModel16 printout. @L479')")
     write(iP,"( 79('/')/)")
   endif !iP>5
   return
End Subroutine DodecModel16
!---------------------------------------------------------------------------7 9
Subroutine SinCos16(AngleD,SinAngle,CosAngle)
!2013.10.22.0940cdt JMS- Using 41! & 40!
!2013.10.14.1150cdt JMS- "Improved" DtoR value confirmed!
!                      - Tested vs. absoft's sinD & cosD in [-720.36,720.36]deg
!                        in .0036 degree steps. Disagreements all below 1.d-18
!2013.10.13.0950cdt JMS- Computes the sine & cosine of an angle
!                      - result based on QPClosure's testing:   2013.10.13.0950
!-----         From the web:   Pi=3.1415926535897932384626433832795028841971...
!   Mine: 29 places:              3.14159265358979323846264338323464000 :Pi
!                     Hence: DtoR=0.0174532925199432957692369076849049
!              Confirmed: SinCos16(30._16,sY,cY)
!              sY= 0.500000000000000000000000000000456   Accurate to 30 places!
!                         SinCos16(60._16,sY,cY)
!              cY= 0.499999999999999999999999999999088   Accurate to 30 places!
!-----
 implicit none                                                       !arguments
  real(16)  ::AngleD   !Input - angle (degrees)
  real(16)  ::SinAngle !Output-   sine(AngleD)
  real(16)  ::CosAngle !      - cosine(AngleD)
  !                                                                  !internals
  real(16)  ::AngleDL
  real(16)  ::DtoR
  real(16)  ::AngleR
  real(16)  ::SumNum
  integer(4)::i,iQ,iQpm,iQuadrant,iQuadrant2
  real(16)  ::Xsq
  real(16)  ::SinAngleL
  !-----
   !-- Internalize the angle:
   AngleDL=AngleD
   !-- Bound to within [-180.,180.] degrees:
   if(AngleDL.gt. 360._16) AngleDL=AngleDL-int( AngleDL/360._16)*360._16
   if(AngleDL.le.-360._16) AngleDL=AngleDL+int(-AngleDL/360._16)*360._16
   if(AngleDL.gt. 180._16) AngleDL=AngleDL-360._16
   if(AngleDL.le.-180._16) AngleDL=AngleDL+360._16
   !-- Determine the quadrant
   iQuadrant=floor(AngleDL/90._16)
   select case (iQuadrant)
    case(-2,2       )  ![-180.,-90.)
     iQuadrant2=-1 ;AngleDL=AngleDL+180._16
    case(-4,-1,0,3,4)  ![-90., +90.]
     iQuadrant2= 0
    case(-3,1       )  !(+90.,+180.]
     iQuadrant2=+1 ;AngleDL=AngleDL-180._16
   end select !iQuadrant
   !-- Convert degrees to radians
  !DtoR  =  0.0174532925199432957692_16   !<- Math Handbook@page 12 to 22places
  !DtoR  =  0.017453292519943295769236907685_16   !JMS previous best
   DtoR  =  0.017453292519943295769236907684905_16  !using Pi{40}/180.

   AngleR=AngleDL*DtoR !Using 28! & 29! ...can the value really be that good?
!                      !Using 36! & 37!        the results were the same.
   Xsq=AngleR*AngleR
   SumNum=0._16    ;iQpm=-1._16  ;              !  41! ~= 3.3452d49
   do i=41,3,-2  ;iQ=i ;iQpm=-iQpm              !  37! ~= 1.3764d43
     SumNum=(SumNum+iQpm)*Xsq/(iQ*(iQ-1._16))   !  29! ~= 8.8417d30
   enddo !i
   SinAngleL=(SumNum+1._16)*AngleR
  !-----
   SumNum=0._16    ;iQpm=-1._16                 !  40! ~= 8.1591d47
   do i=40,2,-2  ;iQ=i ;iQpm=-iQpm              !  36! ~= 3.7199d41
     SumNum=(SumNum+iQpm)*Xsq/(iQ*(iQ-1._16))   !  28! ~= 3.0488d29
   enddo !i
   CosAngle=SumNum+1._16
  !-----
   SinAngle=SinAngleL
   if(iQuadrant2.ne.0) then; Sinangle=-sinangle; CosAngle=-CosAngle ;endif
   return
End Subroutine SinCos16
!---------------------------------------------------------------------------7 9
Subroutine ASinD16(Angle,SineAngle,iP)
!2013.10.20.0810cdt JMS- Computes the asind(SineAngle) to ~31 signif. places.
 use Dodecio, only:sqrj
 implicit none                                                       !arguments
  real(16)  ::Angle    !Output- angle                                    degrees
  real(16)  ::SineAngle!Input - sine of the angle                          "
  integer(4)::iP
  !                                                                  !internals
  real(16)  ::Atemp,sinAtemp,sA,cA
  integer(4)::i,iRev,iSignval
  real(8)   ::DtoR=0.017453292519943295769236907684905_16
  !-----
   sinAtemp=SineAngle    ;iSignVal= 1
   if(sinAtemp.lt.0._16)  iSignVal=-1
   sinAtemp=SineAngle   *iSignVal   ;iRev=0 
   if(sinAtemp.gt..7071_16) then    ;iRev=1 
     sinAtemp=sqrj(1._16-sinAtemp*sinAtemp)
   endif !
   Atemp=dasind(sinAtemp) +.001d0     !getting into the ballpark         
   iRev=0 
      i=0
  !do i=1,10
10   i=i+1
     ATemp=Atemp-(sA-sinAtemp)/DtoR   
     call SinCos16(ATemp,sA,cA)
     if(iP.gt.5) write(iP,"(i2,2f40.33       )")  i,ATemp,(sA-sinAtemp)/DtoR

     if(i.gt.200) then                                              !2024.04.21
       write( 6,"(i3,2f40.33       )")  i,ATemp,(sA-sinAtemp)/DtoR
       pause &
      'pause: ASinD16 @L534 didn`t converge <1.d-30 error in 200 iterations'
       goto 20
     endif!(i.gt.200) 

     if(abs(sA-sinAtemp).gt.1.d-30) goto 10
    !if(abs(sA-sinAtemp).gt.1.d-33) goto 10
  !enddo!i
20 continue
   select case(iRev)
    case(0) ;Angle=        ATemp *iSignVal
    case(1) ;Angle=(90._16-ATemp)*iSignVal
   end select !iRev 
   return
  !Test case:
  !     V1Ang  = -10.812316963571706291284944774472 - previous result
  !sind(V1Ang) =  -0.187592474085079899860139346908 - theoretical
  !     V1Ang  = -10.812316963571706291284944774474 - this algorithm
End Subroutine ASinD16
!---------------------------------------------------------------------------7 9
Subroutine ACosD16(Angle,CosineAngle,iP)
!2013.10.20.0810cdt JMS- Computes the acosd(coainwAngle) to ~31 signif. places.
!                          an Angle in the interval: [0.,+180.] is returned.
 implicit none                                                       !arguments
  real(16)  ::Angle    !Output- angle                                   degrees
  real(16)  ::CosineAngle!Input - cosine of the angle                        "
  integer(4)::iP
  !                                                                  !internals
  real(16)  ::Atemp,cosAtemp
  integer(4)::iSignval
  !-----
   cosAtemp=CosineAngle ;iSignval= 1 
   if(cosAtemp.lt.0.)    iSignval=-1
   call ASinD16(      Atemp,cosAtemp*iSignVal,iP)
                      Atemp= 90._16-Atemp
   if(iSignval.eq.-1) Atemp=180._16-Atemp
                Angle=Atemp
   if(iP.gt.5) write(iP,"(    f40.33)") Angle
   return
  !Test case:
  !     V16Ang = -52.622631859350304357142861505007 - previous result
  !cosd(V16Ang)=   0.607061998206686223095391581420 - theoretical
  !     V16Ang =  52.622631859350304357142861505100 - this algorithm
  !cosd(V16Ang)=   0.607061998206686223095391581419   ... and its cosine
End Subroutine AcosD16
!---------------------------------------------------------------------------7 9
Subroutine hPointPolar16(Ph,aRpyh,PDmag)
!2013.10.17.2210cdt JMS- Quad precision version 
!2013.10.01.2155cdt JMS- Expresses a 3D point in FS ~polar coordinates.
 !-----
 Implicit none                                                       !arguments
  real(16)  ::Ph(4)    !Input - Point#1
  real(16)  ::aRpyh(4) !Output- att.(Roll,Pitch,Yaw)-                  [degrees]
  real(16)  ::PDmag    !      - Distance to point
  !---                                                               !internals
  real(16)  ::Ptemph(4)
  !-----
   Ptemph=Ph
   call hVnorm16(Ph,Ptemph,PDmag)    
   aRpyh=(/0._16,0._16,0._16,1._16/)

   aRpyh(2)=-dasind(Ptemph(3))          !getting into the ballpar            !?
  !call ASinD16(aRpyh(2),Ptemph(3),0)                               !2025.04.20
  !aRpyh(2) =aRpyh(2)

   if((Ptemph(1).ne.0._16).or.(Ptemph(2).ne.0._16))  &
   aRpyh(3)=datan2d(Ptemph(2),Ptemph(1))
   return
End Subroutine hPointPolar16
!---------------------------------------------------------------------------7 9
Subroutine hVnorm16(Vinh,Vouth,Vmag)
!2013.10.17.2210cdt JMS- Quad precision version 
!2013.10.03.2315cdt JMS- Normalizes a homogeneous vector: Vouth=norm(Vinh)
 use Dodecio, only:sqrj
 implicit none                                                        !arguments
  real(16)  ::Vinh(4)    !Vector- input- #1 
  real(16)  ::Vouth(4)   !             - #2
  real(16)  ::Vmag       !      - Dot product = V1h .dot. V2h
  !---                                                                !internals
  real(16)  ::Temph(4)   !Isolates DotP from V1h() & V2h().
  real(16)  ::TempMag
  !-----
                            Temph=Vinh
  if(Vinh(4).ne.0._16) Temph=Temph/Vinh(4)
  TempMag=Temph(1)*Temph(1)+Temph(2)*Temph(2)+Temph(3)*Temph(3)
  if(TempMag.gt.0._16) then
    Vouth(1:3)=Temph(1:3)/sqrj(TempMag) ;Vouth(4)=1._16 ;Vmag=sqrj(TempMag) 
   else
    Vouth=0._16                 ;Vmag=0._16
  endif !TempMag>0
  return
End Subroutine hVnorm16
!---------------------------------------------------------------------------7 9
Subroutine hVecMirror2D16(Vinh,Vmirrorh,Vouth)
!2013.10.09.0800cdt JMS- Mirrors one vector across another
!                          without using trig functions.
 implicit none                                                       !arguments
  real(16)  ::Vouth(4)   !Output- vector  (X,Y,Z)
  real(16)  ::Vmirrorh(4)!Input - mirror vector
  real(16)  ::Vinh(4)    !      - vector
  !                                                                  !internals
  real(16)  ::DotIn
  real(16)  ::VmirrorMagSq
  real(16)  ::ScaleFactor
  real(16)  ::vTemph(4)
  !-----
   DotIn=            Vinh(1)*Vmirrorh(1)  &
                    +Vinh(2)*Vmirrorh(2)  &
                    +Vinh(3)*Vmirrorh(3)
   VmirrorMagSq= Vmirrorh(1)*Vmirrorh(1)  &
                +Vmirrorh(2)*Vmirrorh(2)  &
                +Vmirrorh(3)*Vmirrorh(3)
                             ScaleFactor=2._16*DotIn/VmirrorMagSq
   vTemph(1:3)=Vmirrorh(1:3)*ScaleFactor-Vinh(1:3)
   vTemph( 4 )=                          Vinh( 4 )
   Vouth      =Vtemph
   return
End Subroutine hVecMirror2D16
!---------------------------------------------------------------------------7 9
Subroutine hVecMirror3D16(Vinh,Vmirror1h,Vmirror2h,Vouth)
!2013.10.14.0930cdt JMS- Remirroring back errors are tiny... ~< 1.d-29 
!2013.10.09.0800cdt JMS- Mirrors one vector across two other vectors
!                          without using trig functions.
 implicit none                                                       !arguments
  real(16)  ::Vouth(4)   !Output- vector  (X,Y,Z)
  real(16)  ::Vmirror1h(4)!Input - mirror vector- #1
  real(16)  ::Vmirror2h(4)!                     - #2
  real(16)  ::Vinh(4)    !      - vector
  !                                                                  !internals
  real(16)  ::Vcrossh(4)
  real(16)  ::Vinh2(4)
  !-----
   call hCross16(Vcrossh,Vmirror1h,Vmirror2h)
   call hVecMirror2D16(Vinh,Vcrossh,Vinh2)
   Vouth(1:3)=-Vinh2(1:3)
   Vouth( 4 )= Vinh2( 4 )
   return
End Subroutine hVecMirror3D16
!---------------------------------------------------------------------------7 9
Subroutine hCross16(Voh,Vi1h,Vi2h)
!2013.10.09.1005cdt JMS- Cross product of V1h & V2h
implicit none                                                        !arguments
  real(16)  ::Vi1h(4)     !Input Vector
  real(16)  ::Vi2h(4)     !Input Vector
  real(16)  ::Voh(4)      !Voh = Vi1h X Vi2h
 !--                                                                 !internals
  real(16)  ::Temph(4)
  !----- 
   Temph(1)=Vi1h(2)*Vi2h(3)-Vi1h(3)*Vi2h(2)
   Temph(2)=Vi1h(3)*Vi2h(1)-Vi1h(1)*Vi2h(3)
   Temph(3)=Vi1h(1)*Vi2h(2)-Vi1h(2)*Vi2h(1)
   Temph(4)=Vi1h(4)*Vi2h(4)
   Voh=Temph
   return
End Subroutine hCross16
!---------------------------------------------------------------------------7 9
Subroutine MantExp16(ValIn16,rMant16,iExp,iP)
!2025.04.27.1110cdt JMS- Mantissa and Exponent of a ValIn16
implicit none                                                        !arguments
  real(16)  ::ValIn16     ! Input- value
  real(16)  ::rMant16     !Mantissa
  integer(4)::iExp        !Exponent 
  integer(4)::iP
 !--                                                                 !internals
  real(16)     ::Valu
  character(60)::c60
  character(45)::c45
  !-----
   Valu = ValIn16
   if(Valu>0._16) then
     Write(c60,"(e45.35,    ' :ValTemp')"       ) Valu
     write(c45,"(a39,'E+00'            )"       ) c60(1:39)  
     read( c45,"(e45.35                )"       ) rMant16
     read(c60(42:45),"(i4              )",err=25) iExp
                                         goto 26  
25   read(c60(43:45),"(i3              )"       ) iExp
26   continue
    else
     rMant16 = 0._16
     iExp    = 0
   endif!(Valu>0._16)

   if(iP>5) then
     Write( iP,"(e45.35  ,'  mantissa and exponent:')") Valu
     write(iP ,"(e45.35  ,' :mantissa')"              ) rMant16
     write(iP,"(41x,SP,i4,' :',13x,      'exponent')" ) iExp 
     write(iP ,"( '      12345678901234567890123456789012345')")
     write(iP ,"( '               1         2         3')")
   endif!(iP>5)
    !-----
   return
End Subroutine MantExp16
!         1         2         3         4         5         6
!123456789012345678901234567890123456789012345678901234567890
!    0.25000000000000000000000000000000000E+01  mantissa and exponent:
!    0.25000000000000000000000000000000000E+00 :mantissa
!                                           +1 :             exponent
!---------------------------------------------------------------------------7 9

Subroutine Sqrt16a(ValIn,ValSqrOut,iMZP,iP) !Rev A - 4 arguments.
!2025.04.27.0830cdt JMS- ValSqrOut = sqrt(ValIn)    APF21.0.2: has errors.
!                                                (Perhaps a hardware issue.)
!2013.10.22.0940cdt JMS- This routine isn't needed. APF9.0's : sqrt() is fine!
!2013.10.20.1655cdt JMS- Iteratively computes a real(16) square root.

!The operating domain of the mantissa algorithm is [.1 -to- 10.]
!Odd-order exponents are made even in Iteration i = 1 by:
!    1. multiplying the mantissa by 10._16, and then
!    2. subtracting one from the exponent.
!    3. Setting Y(1) = sqrt(10._16)  ... which brings
!         the iterative domain into [.1 -to- 1.] which is well-behaved.
!The even-order exponent is divided by 2.

implicit none                                                        !arguments
  real(16)  ::ValIn       ! Input- value
  real(16)  ::ValSqrOut   ! = sqrt(abs(ValIn))
  integer(4)::iMZP        !"Minus,Zero, or Plus = {-1,0,+1} =-1:imaginary
!                         ! & ValIn=(ValSqrOut*ValSqrOut)*iMZP
  integer(4)::iP
 !--                                                                 !internals
  real(16)  ::Valu 
  real(16)  ::rMantissa
  integer(4)::iExponent


  real(16)  ::VL,VsqrL,VerrL
  integer(4)::i

  real(16)  ::X(0:9)
  real(16)  ::Y(0:9),Yout
  integer(4)::Init
  character(45)::c45
  real(16)  ::Sqrt10
  !-----
   if(iP>5) write(iP,"(/'Sqrt16(): Valin =')")   !/2x,e43.35)") VaLIn 
   if(     ValIn< 0._16) then; Valu      = -Valin; iMZP = -1
    elseif(ValIn==0._16) then; ValSqrOut =  0._16; iMZP =  0     ;return
    else                     ; Valu      =  ValIn; iMZP = +1
   endif! (ValIn< 0._16)
   call MantExp16(Valu,rMantissa,iExponent,iP)
  !--Mantissa computation:
                   X(0) = rMantissa
   if(iP.gt.5) write(iP,"( '           iMZP =',i3       )") iMZP
   if(iP.gt.5) write(iP,"( 'Mantissa:')")
   Yout = 1._16
   do i = 1, 9;    X(i) = X(i-1)                                   !key ops !
         Y(i) = (((X(i)- 5._16)* X(i)+15._16)* X(i)+ 5._16)/16._16 ! "   "  !
     if(i==1) then  !Accomodatint odd exponents:
       if(mod(iExponent,2).ne.0) then
         Sqrt10 = 3.1622776601683793319988935444327_16 
         Y(i) = Sqrt10
         iExponent = iExponent-1
                   X(0) = X(0)*Sqrt10*Sqrt10 
       endif!(mod(iExponent,2).ne.0)
     endif!(i==1)
     if(iP>5) write(iP,"(i3,2f36.32, ' :i X(',i1,') -> Y(',i1,')')") &
                  i,X(i-1),Y(i),i-1,i
                    X(i) = X(i-1)/(Y(i)*Y(i))
    Yout = Yout*Y(i)
    if(X(0)-Yout*Yout == 0._16) exit
    if(Y(i)           == 1._16) exit
   enddo!i
   if(iP>5) then
               write(iP, &
      "(7x,'Computed Mantissa X value:',10x,'Computed Mantissa Y value:')")
               write(iP,"( i3,2f36.32, ' :i X  -> Y ')")-1,X(0),Yout
               write(iP,"( 3x, f36.32,37x,':  X(0)')") X(0)
               write(iP ,"( '       12345678901234567890123456789012')")
               write(iP ,"( '                1         2         3')")
               write(iP,"( 9x,e30.20,7x, ':Mantissa internal X value discrepancy')") &
                           X(0)-Yout*Yout
   endif!(iP>5)
  !--Exponent computation:
   iExponent = iExponent/2
  !--Conversion to text:
   if(abs(iExponent) <100) write(c45,"(f40.35,'E',SP,i3.2)") Yout,iExponent
   if(abs(iExponent)>=100) write(c45,"(f40.35,    SP,i4.3)") Yout,iExponent
  !--Conversion to real(16) result:
   read(c45,"(e45.35)") ValSqrOut
   if(iP>5) write(iP,"(e45.35,' :ValSqrOut')") ValSqrOut          
   if(iP>5) write(iP,"(i5,40x,' :iMZP  {Minus,Zero, or Plus}')") iMZP 
   if(iMZP==-1) write(iP, &
      "('    ***The square root is an imaginary number.***')")         
   call SaveOutFile
                                                                         return  
End Subroutine Sqrt16a

!Polynomial  coefficients C :  3rd order             ...these are exact!
!1           0.31250000000000000000000000000000 *x^0
!2           0.93750000000000000000000000000000 *x^1
!3          -0.31250000000000000000000000000000 *x^2
!4           0.06250000000000000000000000000000 *x^3
!Y(0) =  (((X(0)-5._16)*x(0)+15._16)*x(0)+5._16)/16._16

! real(16)  ::Xval(-9:18)  & !These are mantissa testcase values of ValIn.
!  =(/  .1_16, .2_16, .3_16, .4_16, .5_16, .6_16, .7_16, .8_16, .9_16 ,1._16  &
!     ,1.1_16,1.2_16,1.3_16,1.4_16,1.5_16,1.6_16,1.7_16,1.8_16,1.9_16         &
!            ,2.0_16,3.0_16,4.0_16,5.0_16,6.0_16,7.0_16,8.0_16,9.0_16, 10._16  /)
!Mantissa processing: 
!  i      Y(1)*Y(1)-X(0)=one-iteration error    X(0) For:
! -9      0.06256016015625000000000000000000     0.1   even order exponents:
! -8      0.03814400000000000000000000000000     0.2
! -7      0.02184347265625000000000000000000     0.3
! -6      0.01152224999999999999999999999999     0.4
! -5      0.00543212890625000000000000000000     0.5
! -4      0.00217600000000000000000000000001     0.6
! -3      0.00067362890625000000000000000001     0.7
! -2      0.00013025000000000000000000000001     0.8
! -1      0.00000797265624999999999999999999     0.9
!  0      0.00000000000000000000000000000000     1.0
!                    1         2         3  
!           12345678901234567890123456789012
!                                                    For:
!  1      0.00000766015624999999999999999999     1.1   odd order exponents:
!  2      0.00012024999999999999999999999999     1.2   forced to even order
!  3      0.00059769140625000000000000000001     1.3   in the 1st iteration
!  4      0.00185600000000000000000000000006     1.4   by Y(1) = sqrt(10)
!  5      0.00445556640625000000000000000000     1.5
!  6      0.00909224999999999999999999999996     1.6
!  7      0.01659128515624999999999999999999     1.7
!  8      0.02790399999999999999999999999996     1.8
!  9      0.04410734765625000000000000000001     1.9
! 10      0.06640625000000000000000000000000     2.0

!           3rd-order correction errors for X(0)>2.  :
! 11      1.00000000000000000000000000000000     3.0 
! 12      5.37890625000000000000000000000000     4.0
! 13     20.00000000000000000000000000000000     5.0
! 14     61.03515625000000000000000000000000     6.0
! 15    162.00000000000000000000000000000030     7.0
! 16    384.53515625000000000000000000000110     8.0
! 17    831.99999999999999999999999999999940     9.0
! 18   1665.87890625000000000000000000000000    10.0
!-----

!Polynomial  coefficients C :  5th order
!1           0.24609375000000000000000000000001 *x^0
!2           1.23046875000000000000000000000000 *x^1
!3          -0.82031250000000000000000000000000 *x^2
!4           0.49218750000000000000000000000000 *x^3
!5          -0.17578125000000000000000000000000 *x^4
!6           0.02734375000000000000000000000000 *x^5

!Y =  ((((( +0.0273437500_16  *X &
!           -0.1757812500_16) *X &
!           +0.4921875000_16) *X &
!           -0.8203125000_16) *X &
!           +1.2304687500_16) *X &
!           +0.2460937500_16)

!  i         Y(1)*Y(1)-X(0)=one-iteration error    X(0)
! -9         0.04255049171447753906250000000000     0.1
!convergence is rapid at 3rd order - hence little improvement:
!-----
!Polynomial  coefficients C :  6th order
!1           0.22558593749999999999999999999999 *x^0
!2           1.35351562500000000000000000000001 *x^1
!3          -1.12792968750000000000000000000000 *x^2
!4           0.90234375000000000000000000000000 *x^3
!5          -0.48339843750000000000000000000000 *x^4
!6           0.15039062500000000000000000000000 *x^5
!7          -0.02050781250000000000000000000000 *x^6

!  i         Y(1)*Y(1)-X(0)=one-iteration error    X(0)
! -9         0.02285984717863880252838134765624     0.1
!convergence is rapid at 3rd order - hence little improvement:
!---------------------------------------------------------------------------7 9

Function Factrl(N) Result(Out)  !<-integer(8)
!2024.07.30.1135cdt JMS- Factorial Computation using a function.
!                          13!=6,227,020,800  > integer(4)
!                          20!=2.432902008e18 < integer(8) is the upper limit.
implicit none                                                        !arguments
  integer(4)::N
  integer(8)::Out
 !--                                                                 !internals
  integer(4)::i
  !-----

   if(N.ge.0) goto 5
     stop 'Factorl: N < 0 '
 5 Out=1_8; if(N.lt.2) return
   do i=2,N; Out=Out*i; end do; return
End Function Factrl
!Used in a print statement:
! integer(8)::Factrl
!   !Test of function Factrl():  2024.07.30:
!    if(iP>5) write(iP,"('Factrl(13) =',i20)") Factrl(13)  <--
!                         Factrl(13) =          6227020800

!2025.04.21.1450cdt JMS- Fortran real(16)'s have ~ a 102-bit mantissa;
!                        hence 2^102 = 5.070602401e30 < 29!
!                        would be the upper limit of mantissa-exact factorials
!                        without more clever coding.
!                        69! = 1.711224524e98
!---------------------------------------------------------------------------7 9

!End Module hVecMath16Mod
!---------------------------------------------------------------------------7 9

