!S3Main.f95          Group ID: #3     Stereo-3D Simulation Environment Vsn:1.00
!2025.05.24.1840cdt- Main program & some utilities.

!            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...
! Program    S3Main
! Subroutine BreakHere(iP,Label)
! Subroutine jPause(Label)
! Subroutine SaveOutFile
! Subroutine M44V4multh(V4outh,Matrix44,V4inh,iP)
! Subroutine M44multh(M44outh,M44Linh,M44Rinh,iP)
! Subroutine h44Fill(Hout,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P)
! Subroutine Invertr8(N,MatrixIn,MatrixOut,iP)
! Subroutine Invertr16(N,A,ValMin,iRank,DetN,iUsed,iP)
! Subroutine Printr16(N,A,Noise,iRank,DetN,iRu,jCu,iP)
! Subroutine Rand16(n,Xrandom,iP)
! Subroutine X12Y12toMB(X1,X2,Y1,Y2,M,B,iP)
! Subroutine Eval(Stage,Descr1,Descr2,Value1,Value2,iP)
! Subroutine DaTime18
! Subroutine RunSeconds
! Subroutine PixelDraw                                               = S18AppF8
! Subroutine nRunSec(nClock,RunSecs,iW)
! Subroutine Jdate20(DaTime)
!---------------------------------------------------------------------------7-9

Program S3Main
!2025.01.03.0900cst JMS-
!--Globals                                             S1ModDef.f95: 2024.03.18
 use TaskDef   ,only: & !Project/Context
     ExeName,ExeVersion,ExeBanner
 use ioDef     ,only: & !Files,Units,TimeStamp,Selfies,Flags
     ExeFileIn,ExeFileOut,EnvNml,EnvIni,LExists &
    ,UserIni,UserNml                            &
    ,Ur,Us,uT,Uread,DaTimeLabel                 &
    ,Pi16          ,ModePi      ,IterPi
 use S2Callback,only:GlutHandoff
 use ViewDef   ,only:FrustNom
 use SimDef    ,only:RunSecs
 use AppsDef   ,only:AppNumber,AppNumberNew

!--End Globals
 implicit none
 !--Arguments
 
 !--Internals
  real(8)   ::RunSecsM
  integer(4)::i
 !real(16)::ValIn,ValSqrOut
  
!S3Main-CallYourApps.f956/ImportYourNML supercedes use of this:
  NAMELIST / Sn3Dnml /  ExeFileIn, ExeFileOut           &
                       ,UserIni  , UserNml              &
                       ,ExeBanner, ExeName, ExeVersion  &
                       ,Ut       , FrustNom

  !--EndDefs-----------------------------------------------
   call nRunSec(0,RunSecsM,6)
!  Modifying key variables at runtime using NameList file `Sn.nml`:
!  Initial defaults are set here:

   if(0>1) then !Bypassed-----------
!  superceded by ImportYourNml() in S3Main-CallYourApps.f95- called below.
   Ut         = 13                  !Output unit#: Set to 6 or 13
   EnvNml     = '.\Sn3D.nml'        !default   
   ExeFileOut = '.\Sn3D-Out.txt'        !   "

   write(Us,"('S3Main: reading ',:,79a1)")(EnvNml(i:i),i=1,len_trim(EnvNml))
   inquire(file=EnvNml,exist=LExists)
   if(LExists)then
     open( unit=Ur, file=EnvNml        , action='read'    &
                  , access='sequential', status='old' ,err=7)
       read(Ur,nml=Sn3Dnml,err=8)     !<- fixed target defined here.
     close(Ur); goto 10
7    pause "Error opening `SnNml.nml` file. Press enter to continue."  ;goto 9
8    pause "Error reading `SnNml.nml` file. Press enter to continue."
9    close(Ur)
    else
     write(Us,*) "Didn't find `SnNml.nml`."
   endif!LExists
10 continue
!  write(Us,"('SnNml.nml: UserNml(2) =',:,79a1)") &
!                     (UserNml(2)(i:i),i=1,len_trim(UserNml(2))) !Diagnostic
   endif!(!(0>1) end Bypassed-----
   
   call ImportYourNml !In S3Main-CallYourApps.f95

   write(6,"('S3Main: ExeFileOut     = ',:,79a1)") &
                     (ExeFileOut(i:i),i=1,len_trim(ExeFileOut))
   if(Ut<10) goto 18
   open(unit=Ut,file=ExeFileOut,action='write',err=19); goto 20
18 Ut=Us; write(Us,"(/'You disabled the output text file by setting Ut<10')")
          write(Us,"('Ut=Us; output text -> DOS screen.')")
          goto 20
19 Ut=Us; write(Us, &
     "(/'Unable to open the output text file (e.g. running from a CD)')")
          write(Us,"('Ut=Us; & output text -> DOS screen.')")
20 call DaTime18
   write(Ut,"(a18,' Begin.',5x,a22,1x,a22)") DaTimeLabel,ExeName,ExeVersion

  !Call/test a non-visual subroutine here:
   if(0>1) then !---------- Bypasses when (0>1); enters using (1>0)
    !call ________ 
     call SaveOutFile    !...this will save the printout.
     stop'S3Main @L1117'
   endif!(0>1) -------- end bypass

  !Perform your application-specific pre-OpenGL initialization(s), if any:
  ! call _________
  !pre-OpenGL initializations done.
   if(AppNumber<1) AppNumber = 1 
   if(AppNumber>4) AppNumber = 1 

   call GlutHandoff  !OpenGL will take control...
End Program S3Main
!---------------------------------------------------------------------------7 9

Subroutine BreakHere(iP,Label)
!2020.05.10.0505cdt JMS- Call this routine from a problem area in your code.
!  call BreakHere(iP,'Label') !////////////////////////////////////////////////
!  call BreakHere(iP)         !////////////////////////////////////////////////
 implicit none
 !--Arguments
  integer(4)                       ::iP    !Write enable>5: write(iP,...)
  character(*),intent(in),optional::Label !Ref.automatic: P.109
 !--Internals
  integer(4)::i,j
  !--EndDefs-----------------------------------------------
   if(iP>5) then
     if(present(Label)) then     ;write(iP,"('BreakHere:',\)")
       j=len_trim(Label); if(j>0) write(iP,"( 80a1:)") (Label(i:i),i=1,j)
     endif!present: Label
   endif!iP>5
   return  !<- Breakpoint. From here, you'll return to the calling location.
End Subroutine BreakHere
!---------------------------------------------------------------------------7 9

Subroutine jPause(Label)
!2022.04.19.1340cdt JMS- Replaces the depricated "Pause".
 implicit none                                                       !arguments
  character(*),intent(in)::Label  !Label to write to unit# iWU, ="":no-op
  !---                                                               !internals
  integer(4)   ::i,j
   !-----                                                              end defs
!  Determine the Label length:
   j = len_trim(Label)
!  Write Label + " Press `Enter` to continue: "
   if(j>0) write(6,"(80a1: )",advance="no") (Label(i:i),i=1,j)
           write(6,"(' Press `Enter` to continue: ')",advance="no")
!  The pause is waiting for a keyboard return:
   read(5,*)   !Keyboard attended, hence pausing
   return
end Subroutine jPause
!---------------------------------------------------------------------------7 9

Subroutine SaveOutFile                                         !2023.00.02.0640
!2023.06.24.1155cdt JMS
!--Globals
 use ioDef,only:ExeFileOut,uT
!--End Globals
 implicit none
 !--Arguments
 !--Internals
  !--EndDefs-----------------------------------------------
   if(Ut<10)                                                             return
   close(Ut)
   open(Ut,file=ExeFileOut,action='write',position='Append')            ;return
End Subroutine SaveOutFile
!---------------------------------------------------------------------------7-9

Subroutine M44V4multh(V4outh,Matrix44,V4inh,iP)
!2020.05.05.1745cdt JMS- V4outh(4) = Matrix44(4,4) * V4inh(4)
 implicit none
 !--Arguments
  real(8)   ::V4outh(  4  )
  real(8)   ::Matrix44(4,4)
  real(8)   ::V4inh(   4  )
  integer(4)::iP            !Write enable>5: write(iP,...)
 !--Internals
  integer(4)::i,j
  real(8)   ::Vtemp(4)
  !--EndDefs-----------------------------------------------
   do i=1,4
     Vtemp(i) =  Matrix44(i,1)*V4inh(1)  &
               + Matrix44(i,2)*V4inh(2)  &
               + Matrix44(i,3)*V4inh(3)  &
               + Matrix44(i,4)*V4inh(4)
   enddo!i
   Vtemp=Vtemp/Vtemp(4)
   if(iP>5) then
     write(iP,"(5x,'V(4)out',27x,'M(4,4)',26x,'V(4)in')")
     write(iP,"(f12.6,'   ',4f12.4,'   ',f12.6)")  &
                Vtemp(1),(Matrix44(1,j),j=1,4),V4inh(1)
     write(iP,"(f12.6,' = ',4f12.4,' * ',f12.6)")  &
                Vtemp(2),(Matrix44(2,j),j=1,4),V4inh(2)
     write(iP,"(f12.6,'   ',4f12.4,'   ',f12.6)")  &
                Vtemp(3),(Matrix44(3,j),j=1,4),V4inh(3)
     write(iP,"(f12.6,'   ',4f12.4,'   ',f12.6)")  &
                Vtemp(4),(Matrix44(4,j),j=1,4),V4inh(4)
   endif !iP>5

   V4outh=Vtemp
   return
End Subroutine M44V4multh
!---------------------------------------------------------------------------7 9

Subroutine M44multh(M44outh,M44Linh,M44Rinh,iP)
!2025.05.02.1255cdt JMS- M44outh(4,4) = M44Linh(4,4) * M44Rinh(4)
 implicit none
 !--Arguments
  real(8)   ::M44outh(4,4)
  real(8)   ::M44Linh(4,4)
  real(8)   ::M44Rinh(4,4)
  integer(4)::iP            !Write enable>5: write(iP,...)
 !--Internals
  integer(4)::i,j
  real(8)   ::M44temp(4,4)
  !--EndDefs-----------------------------------------------
   do i=1,4
     do j=1,4
       M44temp(i,j) =  M44Linh(i,1)*M44Rinh(1,j)  &
                     + M44Linh(i,2)*M44Rinh(2,j)  &
                     + M44Linh(i,3)*M44Rinh(3,j)  &
                     + M44Linh(i,4)*M44Rinh(4,j)
     enddo!j
   enddo!i
   if(iP>5) then
     write(iP,"(4x,'M44outh(4,4)',37x,'= M44Linh(4,4)',33x,'+ M44Rinh(4,4')")
     do i=1,4
       write(iP,"(4f12.6,' ',4f12.4,' ',4f12.6)")  &
          (M44temp(i,j),j=1,4),(M44Linh(i,j),j=1,4),(M44Rinh(i,j),j=1,4)
     enddo!i
   endif !iP>5

   M44outh = M44temp
   return
End Subroutine M44multh
!---------------------------------------------------------------------------7 9

Subroutine h44Fill(Hout,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P)
!2020.04.08.2055cdt JMS- homogeneous(4,4)= A row-order set of 16 coefficients.
 implicit none        !- fills a matrix which has column-order coefficients.
 !--Arguments
  real(8)   ::Hout(4,4)            !The 4x4 matrix to fill
  real(8)   ::A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P !The coefficients of Hout
 !integer(4)::iP                   ! >5: Printing enable, use: write(iP,...)
 !--Internals
  !--EndDefs-----------------------------------------------
   Hout(1,1:4)=(/ A,B,C,D /)
   Hout(2,1:4)=(/ E,F,G,H /)
   Hout(3,1:4)=(/ I,J,K,L /)
   Hout(4,1:4)=(/ M,N,O,P /)                                            ;return
End Subroutine h44Fill
!---------------------------------------------------------------------------7-9
Subroutine Invertr8(N,MatrixIn,MatrixOut,iP)                    !2020.02.24 JMS
!2020.05.01.0550cdt JMS- A(0:15,0:30) <-correction
 !MatrixOut(N,N) = Inverse[MatrixIn(N,N)] by sequential column processing.
 !                  Homogeneous matrices invert and multiply 
 !                        just like other matrices... 
 !                         an amazingly useful fact!
 !"Matrix inversion" is hyperspace's generalization of division.
 !                     If: [Y]=        [A]*[X]
 !                   Then: [X]=    [Y]/[A]     doesn't exist,  -but-
 !                         [X]= Inverse[A]*[Y] is defined      -but-
 !  [A] has to be "full rank" to invert, which is like being non-zero in 1D.
 !Advanced topics: 
 !     >Inverters more powerful than "Invertr8" provide partial inverses
 !         when matrices are less than "full rank" but not all zero's.
 !     >"Pseudoinverters" provide "least-squares best fits" of arbitrarily
 !         large datasets, and very few lines of code can accomplish that.
 !-----
 implicit none                                                       !arguments
  integer(4)::N
  real(8)   ::MatrixIn( N,N)
  real(8)   ::MatrixOut(N,N)
  integer(4)::iP            !Write enable>5: write(iP,...)
 !---                                                                !internals
  real(8)   ::A(0:15,0:30)  !This intermediate array supports 1<= N <=15
                           !The 0th row is swap space.
  integer(4)::i,iUse,j,Jouter,k
  real(8)   ::Smallest,Largest,Temp,Det
 !character::Msg*79
 !character::Msg*50
  !-----
  !if(N< 1) then; Msg="sub:Invertr8: N< 1"; call AbortNow(Msg); endif !n< 1
  !if(N>15) then; Msg="sub:Invertr8: N>15"; call AbortNow(Msg); endif !n>15
   if(N< 1) stop "sub:Invertr8: N< 1"                                 !n< 1
   if(N>15) stop "sub:Invertr8: N>15"                                 !n>15
             A         = 0.d0          !Clear the working arrays.
             A(1:N,1:N)= MatrixIn      !Import the input matrix.
   do i=1,N; A(i  ,i+N)= 1.d0;enddo!i  !Appends an identity matrix

  !Find the largest abs value:
   Largest=0.d0
   do i=1,N;do j=1,N; if(dabs(A(i,j))>Largest) Largest=dabs(A(i,j));enddo;enddo
   if(Largest==0.d0)then 
     !Msg="sub:InvertMM: input was all zero's."; call AbortNow(Msg)
     !Msg="sub:InvertMM: input was all zero's."; stop Msg
     stop "sub:InvertMM: input was all zero's."
   endif !Largest==0.
  !Set the noise floor at a billionth of the largest matrix coefficient:
   Smallest=Largest/1.d9
   Det=1.d0

  !Solve
                                                       iUse=-1
   do Jouter=0,N   !Outer loop:
     if(Jouter==0) goto 100

     !Find the largest remaining coefficient.
                            Largest=0.d0             ; iUse= 0
     do i=Jouter,N
       if(dabs(A(i,Jouter))>Largest) then
                            Largest=dabs(A(i,Jouter)); iUse= i
       endif !A(i,j)>Largest
     enddo!i
     if(iUse==0) goto 90 !MatrixIn is less than full rank.

     Largest=A(iUse,Jouter) !The signed value.
     Det=Det*Largest
     !Divide the "used row" by the "Largest" value:
     A(iUse,1:2*N)=A(iUse,1:2*N)/Largest

    !If necessary - swap rows:
     if(iUse/=Jouter) then
       A(   0  ,0:2*N)=A( iUse ,0:2*N)
       A( iUse ,0:2*N)=A(Jouter,0:2*N)
       A(Jouter,0:2*N)=A(   0  ,0:2*N)
     endif !iUse/=Jouter

     if(iP>5) then
       write(iP, &
       "( 'Process column',i3,',  largest was',F12.6,' (',i2,',',i2,')')")  &
                   Jouter    ,  Largest        ,   iUse  ,Jouter
       write(iP,"(21x,'Det=',f24.12)") Det
       do i=1,N
         write(iP,"(30f12.6)") (A(i,j),j=1,2*N)
       enddo!i
     endif !iP>5

    !Reduce the Jouter column to an identity.
     do i=1,N; if(i.eq.Jouter) cycle
                                 Temp=A( i    ,Jouter)
       do j=1,2*N; A(i,j)=A(i,j)-Temp*A(Jouter, j    ); enddo!j
     enddo!i
     goto 100 !Inversion proceeding normally...

90   continue !Full inversion failed...
    !Report & exit...

100  continue
     if(iP>5) then  !Report progress:
       if(Jouter==0) &
         write(iP,"(/'sub:Invertr8, N=',i2,'  initial setup:',45('-'))") N
       if(Jouter> 0) write(iP,"( 'Iteration:',i2,' done:')") Jouter   
       do i=1,N
         write(iP,"(30f12.6)") (A(i,j),j=1,2*N)
       enddo!i

       if(Jouter==N) then !Final pass: evaluate the quality of the result.
         write(iP,"( 'Residual Errors of [MatrixOut]*[MatrixIn]-[I]:',\)")
         write(iP,"( '    (All zero`s would be ideal.)')")  
         do i=1,N
           do j=1,N; Temp=0.d0; if(i==j) Temp=Temp-1.d0
             do k=1,N; Temp=Temp+A(i,N+k)*MatrixIn(k,j); enddo!k
             write(iP,"(d15.6,\)") Temp
           enddo!j
             write(iP,*)
         enddo!i
         write(iP,"('sub:Invertr8 done.')")
       endif !Jouter==N
     endif !iP>5

     if(iUse==0) then
      !Msg="sub:Invertr8: Failed. MatrixIn wasn't invertable"
      !call AbortNow(Msg)
       stop "sub:Invertr8: Failed. MatrixIn wasn't invertable"
     endif !n< 1
    !Inversion proceeding
   enddo!Jouter

   MatrixOut=A(1:N,N+1:2*N)
   return
End Subroutine Invertr8
!---------------------------------------------------------------------------7 9

Subroutine Invertr16(N,A,ValMin,iRank,DetN,iUsed,iP)
!2010.01.06.1205cst JMS- Traveler2/Athlon64/WinXPPro-32/APF9.0-32
!2006.05.25.0715cdt JMS- An instantiation of the "Setterholm Matrix Inverter"
!                          Expects,uses, & returns real(16)'s (~QUAD precision).
!---                                                                 No globals
 implicit none                                                       !arguments
  integer(4):: N
  real(16)  :: A(N,N)
  real(16)  :: ValMin
  integer(4):: iRank
  real(16)  :: DetN    !Abs(determinant product) for each iteration
  integer(4):: iUsed(N)
  integer(4):: iP
  !---                                                               !internals
  integer(4),allocatable::iRu(:)
  integer(4),allocatable::jCu(:)
  real(16)  ,allocatable::Temp(:)
  integer(4)::i,iu,iu2,j,ju,ju2,L
  real(16)  ::Amult,ValMax
  integer(8)::iBinSum
  integer(4)::iAlloc
   !-----                                                            !end defs
   if(iP.gt.5) write(iP,"('Overwriting Matrix Inversion:')")
   DetN=0._16 ;iUsed(1:N)=0 ;iRank=0
   !---
   allocate(iRu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: iRu(N) allocation error. Halt.'
   iRu=0
   allocate(jCu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: jCu(N) allocation error. Halt.'
   jCu=0
   allocate(Temp(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'Invert: Temp(N) allocation error. Halt.'
   Temp=0._16
   !---   
     Temp=0._16; iRu   =0 ;jCu   =0 ;ValMin=0._16
   do i=1,N    ;iRu(i)=i ;jCu(i)=i                ;end do
           if(iP.gt.5) call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
!  Inversion Loop:
   do L=1,N ;iU=0     ;jU=0     ;ValMax =0._16
     do i=1,N  ;if(iRu(i).lt.0) cycle
       do j=1,N ;if(jCu(j).lt.0) cycle
         if(ValMax.ge.abs(A(i,j))) cycle
           ValMax=abs(A(i,j)) ;iu=i ;ju=j
       enddo!j
     enddo!i
     if(L==1) ValMin=ValMax*1.d-24  !Establishes the "noise floor". Was 1.d-20
     if(L==1) DetN=ValMax
     if(iU.eq.0) Exit
     if(ValMax.lt.ValMin) Exit
     iRank=iRank+1
     if(L>1) DetN=DetN*ValMax
     ValMax=A(iu,ju) ;iu2=jCu(ju) ;ju2=iRu(iu)
    !if(iP.gt.5) write(iP,"(3i3,f15.6)")  L,iu,ju,ValMax
     if(iP.gt.5) write(iP,"(1x,2i4,f15.9,'   :iu ju pivot, & value')") &
                                              iu,ju,        ValMax
     Temp(1:N)=A(1:N,ju)                            ;i=iRu(iu)
               A(1:N,ju)=A(1:N,ju2)                 ;  iRu(iu)=iRu(iu2)
                         A(1:N,ju2)=Temp(1:N)       ;          iRu(iu2)=-abs(i)
     Temp(1:N)=A(iu,1:N)                            ;j=jCu(ju)
               A(iu,1:N)=A(iu2,1:N)                 ;  jCu(ju)=jCu(ju2)
                         A(iu2,1:N)=Temp(1:N)/ValMax ;         jCu(ju2)=-abs(j)
     do i=1,N ;if(i.eq.iu2) cycle ;Amult=A(i,ju2)
       do j=1,N ;A(i,j)=A(i,j)-A(iu2,j)*Amult ;enddo!j
       A(i,ju2)=-Amult/ValMax
     enddo; A(iu2,ju2)=1._16/ValMax
           if(iP.gt.5) call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
   enddo!L  Inversion loop- done.

   iBinSum=0
   do i=1,n
     if((iRu(i).lt.0).and.(i.lt.63)) iBinsum=iBinSum+2**(i-1)
      if(iRu(i).lt.0) iUsed(i)=1
!     Zero linearly-dependent rows and columns, if any:
      if(iRu(i).gt.0) A(i,1:N)=0._16
      if(jCu(i).gt.0) A(1:N,i)=0._16
   enddo!i
   if((iP >5).and.(iu == 0)) then
     write(iP,"('Linear dependency dealt with. Partial inverse results:')")
     write(iP,"(b63.2,' -Kp unused flag')") iBinSum !For the first 62 Kp's
                       call Printr16(N,A,ValMin,iRank,DetN,iRu,jCu,iP) !*******
   endif!iP>5 & iu==0
   if(iP.gt.5) write(iP,"('Overwriting Inverter- done.')")
   deallocate(iRu)
   deallocate(jCu)
   deallocate(Temp)
   return
End Subroutine Invertr16
!---------------------------------------------------------------------------7-9

Subroutine Printr16(N,A,Noise,iRank,DetN,iRu,jCu,iP)
!2018.09.04.1340cdt JMS- Prints Invert's progress.
!                          Expects,uses, & returns real(16)'s (~QUAD precision).
 !---                                                                No globals
 implicit none                                                       !arguments
  integer(4):: N
  real(16)  :: A(N,N)
  real(16)  :: Noise
  integer(4):: iRank
  real(16)  :: DetN
  integer(4):: iRu(N)
  integer(4):: jCu(N)
  integer(4):: iP
  !---                                                               !internals
  integer(4)::i,j
   !-----                                                            !end defs
                                                              if(iP < 6) return
   write(iP,"('Inverter: @ Rank =',i5,8x,'abs(Det)=',e39.30/2x,\)") iRank,DetN
   do j=1,N;   write(iP,"(i18,\)") jCu(j) ; enddo; write(iP,"('')")
   do i=1,N;   write(iP,"(i5,' ' ,\)") iRu(i)
     do j=1,N; write(iP,"(sp,f18.9,\)") A(i,j); enddo; write(iP,"('')")
   enddo!i
   write(iP,"('')");                                                     return
End Subroutine Printr16
!---------------------------------------------------------------------------7-9

Subroutine Rand16(n,Xrandom,iP)                                !2015.01.31.1645
!2015.01.31.1645cst JMS- Random number generator
!                      - Set Xrandom=0._16 prior to first call.
!-- Use:-------------------------- 
! integer*4::nRandom=12
! real*16  ::Xrandom(12)
!Xrandom(1)=-1408701324  -or-
!Xrandom(1)=0._16
!do i=1,10
!  call Rand16(nRandom,Xrandom,iP)
!enddo!i--------------------------
implicit none
 !--Arguments
  integer(4)::n
  real(16)  ::Xrandom(n)
  integer(4)::iP
 !--Internals
  integer(4)::Init
  integer(4)::nRand    !Number of integers in the seed
  integer(4)::iSeed(4)
  integer(4) time
  !--EndDefs-----------------------------------------------
  
 !-- Random number generation initialization:
   if(Init.eq.0) then
     call random_seed
     call random_seed(size=nRand)           !=  1
     call random_seed(get=iSeed(1:nRand))   != -1408701324 (default value)
     if(iP.gt.5) write(iP,"('Rand16: nRand=',i2,' iSeed=',10(i14:))")  &
                                     nRand     ,  iSeed(1:nRand)
       iSeed=0
    !-- Initialize randomness...
     if(Xrandom(1).eq.0._16) then
      !-- ... using the computer's seconds clock when Xrandom(1)=0.:
       iSeed(1)=time()                 !e.g.:=  1422733593 @ 2015.01.31.1346.33
                                       !     =  1422733725 @ 2015.01.31.1348.46
                        !Returns the seconds since 00:00:00 GMT January 1. 1970
      else 
      !-- ... by you setting the integer value(s) of Xrandom(1:nRand):
       iSeed(1:nRand)=Xrandom(1:nRand)+sign(.1_16,Xrandom(1:nRand))
     endif !

     call random_seed(put=iSeed(1:nRand))
     if(iP.gt.5) write(iP,"(5x,'Value used: iSeed=',10(i14:))") iSeed(1:nRand)
     Init=1
   endif !Xrandom(1)=0.
  
   call random_number(Xrandom)
!  For iSeed(1)=-1408701324 (default value)
!  & n=2:
!  1  0.1390874385833740234375000  0.6339781284332275390625000
!  2  0.2369287014007568359375000  0.4304550290107727050781250
!  & n=1:                            123456789012345678901^granular beyond ~21
!  1  0.1390874385833740234375000
!  2  0.6339781284332275390625000  < ^ :Exactly the same sequence.
!  3  0.2369287014007568359375000
!  4  0.4304550290107727050781250
   if(iP.gt.5) then
     if(Init.le.999) then; write(iP,"('Random# set:',i3 ,\)") Init
                     else; write(iP,"('R#:'         ,i12,\)") Init
     endif!Init<=999
     if(n.lt.3) then; write(iP,"(     2(f26.22:)   )") Xrandom(1:n)
                else; write(iP,"(667( /3(f26.22:) ) )") Xrandom(1:n)
!Rand16: nRand= 1 iSeed=   -1408701324
!     Value used: iSeed=    1422743635
!Random# set:  1
!  0.6806479096412658691406  0.6724007725715637207031  0.0493627786636352539063
!  0.4378559589385986328125  0.2765567302703857421875  0.4999549388885498046875
!  0.3916140198707580566406  0.3895425200462341308594  0.3124939203262329101563
!  0.6459279060363769531250  0.5968738198280334472656  0.4783441424369812011719
!Random# set:  2
!  0.7516865134239196777344  0.2383043169975280761719  0.4446819424629211425781
!  0.7403323650360107421875  0.0192565321922302246094  0.0300652384757995605469
!  0.5761010050773620605469  0.7211880683898925781250  0.7397087216377258300781
!  0.9435601830482482910156  0.7619655728340148925781  0.2017260193824768066406
     endif !n<3
   endif !iP>5
   Init=Init+1
   return
End Subroutine Rand16
!---------------------------------------------------------------------------7-9

Subroutine X12Y12toMB(X1,X2,Y1,Y2,M,B,iP)
!2020.05.13.0840cdt JMS- Slope & intercept defined by a line's endpoints.
 implicit none
 !--Arguments
  real(8)   ::X1,X2,Y1,Y2 !Line endpoints
  real(8)   ::M,B         !Slope & Intercept   Y=M*X+B
  integer(4):: iP         !Write enable>5: write(iP,...)
 !--Internals
  !--EndDefs-----------------------------------------------
   M=(Y2-Y1)/(X2-x1); B=Y1-M*X1                                ;if(iP<6) return
  !--Print results to unit iP
   write(iP,"(/'X12Y12toMB: X2=',f15.6,' Y2=',f15.6)") X2,Y2
   write(iP,"( '            X1=',f15.6,' Y1=',f15.6)") X1,Y1
   write(iP,"( '    Deltas:    ',f15.6,'    ',f15.6)") (X2-X1),(Y2-Y1)
   write(iP,"( '     Hence: M =',f15.6,' B =',f15.6)") M,B                ;return
End Subroutine X12Y12toMB
!---------------------------------------------------------------------------7-9

Subroutine Eval(Stage,Descr1,Descr2,Value1,Value2,iP)
!2020.05.07.1345cdt JMS- Incremental numerical derivation confirmation.
 implicit none
 !--Arguments
  character(*),intent(in)::Stage,Descr1,Descr2            !Ref.automatic: P.109
  real(8)   ::Value1,Value2
  integer(4)::iP                   ! >5: Printing enable, use: write(iP,...)
 !--Internals
  integer(4)::i
  !--EndDefs-----------------------------------------------
   if(iP<6)                                                              return
   write(iP,"(/'Eval Stage: ',60a1)") (Stage(i:i),i=1,len_trim(Stage))
   write(iP,"(f19.9,' = ',60a1: )") Value1,(Descr1(i:i),i=1,len_trim(Descr1))
   write(iP,"(f19.9,' = ',60a1: )") Value2,(Descr2(i:i),i=1,len_trim(Descr2))
   return
End Subroutine Eval !Usage example:
  !--        ,"                                             60 char limit-> " &
!  call Eval( "#1->#5" &
!            ,"-( (b+1)*(2*d+m)+4*d-(b-1)*(m-2*d) )" &
!            ,"-(4*b*d+4*d+2*m)"                     &
!            , -( (b+1)*(2*d+m)+4*d-(b-1)*(m-2*d) )  &
!            , -(4*b*d+4*d+2*m)                    ,iP)

!Eval Stage: #1->#5
!     -236.210494995 = -( (b+1)*(2*d+m)+4*d-(b-1)*(m-2*d) )
!     -236.210494995 = -(4*b*d+4*d+2*m)!Eval Stage: #1->#5
!---------------------------------------------------------------------------7-9

Subroutine DaTime18
!2020.04.05.1510cdt JMS- Present date & time string -> DaTimeLabel*18
!--Globals                                             S1ModDef.f95: 2024.03.18
 use ioDef,only:DaTimeLabel !Input/Output & Flags
!--End Globals
 implicit none
 !--Internals
  integer(4)::i,iDMY(3),iHMS(3)
  !--EndDefs-----------------------------------------------

   call iDate(iDMY)
   call iTime(iHMS)
   write(DaTimeLabel,"( i4,'.',i2.2,'.',i2.2,'.',i2.2,i2.2,'.',i2.2)") &
                      (iDMY(i),i=3,1,-1),(iHMS(i),i=1,3)                ;return
End Subroutine DaTime18
!---------------------------------------------------------------------------7 9

Subroutine RunSeconds
!2020.06.16.0640cdt JMS- Run seconds.                             Not verified.
!--Globals
 use SimDef,only:RunSecs   !Input/Output & Flags
!--End Globals
 implicit none
 !--Internals
  integer(4)::Init=0
  integer(4)::Vo(8),V(8)
  integer(4)::V3prev,iDays=0
! integer(4)::DayWrap=0
  !--EndDefs-----------------------------------------------
   if(Init==0) then
     call date_and_Time(VALUES=Vo)
     V=Vo
     Init=1
   endif!Init=0
   V3Prev=V(3)
   call date_and_Time(VALUES=V)
   if(V(3)/=V3Prev) iDays=iDays+1
   RunSecs = ((iDays*24+V(5)-Vo(5))*  60.d0 & !Hours
                       +V(6)-Vo(6))*  60.d0 & !Minutes
                       +V(7)-Vo(7)          & !Seconds
                      +(V(8)-Vo(8))/1000.d0   !Millisec
   return
End Subroutine RunSeconds
!---------------------------------------------------------------------------7 9

!This is example app F8:

Subroutine PixelDraw                                    !2020.02.04.0740cst JMS
!Draws magnifiable pixels directly onto the OpenGL screen.
!  Support seeing pixel-level algorithms running live step-by-step.
!use OglCanMod !<-  S.*, OpenGL, & the environment's subroutines     ! globals
 !-----

 use OpenGLRec              !Ref: OpenGL GL/GLU/GLUT docs
 use ioDef    ,only: Up     !Files,Units,TimeStamp,Selfies,Flags
 use ScreenDef, only: xWindowFull,yWindowFull,ColorRGBA
 use AppsDef  ,only: AppBanner,AppInFile
 use ViewDef  ,only: ThreePhase !View bounds->Unit Cube ->[extent]Destination
!use OglCanDef  !, only: S
!use OglCanMod, only:Colors3D,PrntOrtho
 implicit none                                                       !arguments
  !--                                                                !internals
  integer(4)::MtxMode(1)
  real(8)   ::Orth6(6)
  integer(4)          ::i
  real( kind=glFloat)::w0,h0, PixMagf
  real( kind=glFloat)::fImage(3)=(/ 0., 1., 1. /)
  real(8)             ::PixMag=32
  integer(4)::iColor
  character(len=80)::PText
   !-----                                                            ! end defs
   AppBanner(8)  ='PixelDraw                     2020.02.04.1800'  !a58
   AppInFile(8)  = char(0)                                         !a35
   if(Up>5)write(Up,"('F8 PixelDraw: ThreePhase=',i2,' Up=',i2)")ThreePhase,Up
   select case(ThreePhase)
    case(1)  !Background number crunching, if any:
    case(2)  !Update variables & 2D screen info: Superceded by split screen(s)
    case(3)  !Draw/redraw app. 2D & 3D graphics:

   call glGetIntegerv(GL_MATRIX_MODE,MtxMode)
   call glMatrixMode(GL_MODELVIEW ) ;call glPushMatrix ;call glLoadIdentity
   call glMatrixMode(GL_PROJECTION) ;call glPushMatrix ;call glLoadIdentity

   call glMatrixMode(GL_PROJECTION) ;call glLoadIdentity
   call glMatrixMode(GL_MODELVIEW)  ;call glLoadIdentity

   iColor=0
   Orth6=(/ 0.d0                ,  xWindowFull/PixMag  &
           ,  yWindowFull/PixMag,     0.d0             &
           ,    -1.d0           ,     1.d0              /) 
   call glOrtho( Orth6(1),Orth6(2),Orth6(3),Orth6(4),Orth6(5),Orth6(6) )

   !Alter screen pixel (w0,h0):
                            PixMagf=PixMag
   call glPixelZoom(PixMagf,PixMagf)
   do i=1,495/idint(PixMag)
     iColor=iColor+1; if(iColor.gt.16) iColor=1
     call Colors3D(iColor); FImage=ColorRGBA(1:3)

                       ! w0=i; h0=i*2.0
                         w0=0; h0=i
     call glRasterPos2f( w0  , h0  )
     call glDrawPixels(1_GLsizei,1_GLsizei,GL_RGB,GL_FLOAT,fImage)

     write(Ptext,"('Color#:',i3)") iColor
                                    call PrntOrtho(iColor+4,50,iColor, 0,PText)

   enddo!i
   call glFlush

16 call glMatrixMode(GL_MODELVIEW ) ;call glPopMatrix
   call glMatrixMode(GL_PROJECTION) ;call glPopMatrix
   call glMatrixMode(MtxMode(1))

   end select!ThreePhase
   return
End Subroutine PixelDraw
!---------------------------------------------------------------------------7 9

Subroutine nRunSec(nClock,RunSecs,iW)
!2023.08.12.1645cdt JMS- Run seconds for (0:30) separate clocks
!--Globals                                      ... 1. millisecond resolution
 use ioDef, only: DaT,DaTo
!--End Globals
 implicit none
 !--Arguments
  integer(4)  ::nClock   !<0:reset&start; >0:read; =0:time since first call
  real(8)     ::RunSecs  !nClock`s elapsed seconds
  integer(4)  ::iW
 !--Internals
  integer(4)  ::Init = 0
  integer(4)  ::nClockL
  character   ::DaTime*20
  !--EndDefs-----------------------------------------------
   call Jdate20(DaTime) !Computes DaT
   nClockL = nClock
   if(Init==0) then
     DaTo(0) = DaT
     Init    = 1
   endif!(Init==0)
   select case(nClock)
    case(   :-31); Stop 'nRunSec(nClock <-30... undefined. Halt@L125'
    case(-30:- 1); nClockL = -nClockL; DaTo(nClockL) = DaT
    case(+31:   ); Stop 'nRunSec(nClock >+30... undefined. Halt@L127'
   end select!(nClock)
   RunSecs =  (DaT%Sec -DaTo(nClockL)%Sec )          &
             +(DaT%JDay-DaTo(nClockL)%JDay)*86400.d0
  !if(iW>5) write(iW,"( 5x,'RunSecs(',i2,')=',f27.3)") nClock,RunSecs
   if(iW>5) write(iW,"(10x,SP,i3,S,f31.3)") nClock,RunSecs
                                                                         return
End Subroutine nRunSec
!---------------------------------------------------------------------------7 9

Subroutine Jdate20(DaTime)                                     !2023.08.20.1210
!2023.08.11.0600cdt JMS- Reads a snapshot of the computers current date & time.
 use ioDef, only:DaT
 implicit none
 !--Arguments
  character::DaTime*20
 !--Internals
  integer(4)::V(8)
  integer(4)::iYear,iMonth,iDay,iHour,iMin,iSec,imSec
 !integer(4)::JulianDay
 !real(8)   ::DaySecond
  !--EndDefs-----------------------------------------------
   call date_and_Time(VALUES=V)  !V(4) =Time offset w.r.t. UTC
   iYear  = V(1); iHour = V(5)
   iMonth = V(2); iMin  = V(6)
   iDay   = V(3); iSec  = V(7)
                  imSec = V(8)     !includes:      seconds
   write(DaTime,"(i4,'.',i2.2,'.',i2.2,'.',i2.2,i2.2,'.',i2.2,' L'    )")  &
                  iYear,iMonth,iDay,iHour,iMin,iSec
   write(DaT%cDaT22, &             !includes: milliseconds
                "(i4,'.',i2.2,'.',i2.2,'.',i2.2,i2.2,'.',i2.2,'.',i3.3)")  &
                  iYear,iMonth,iDay,iHour,iMin,iSec,imSec
  !JulianDay
   DaT%jDay =  367*  iyear                                       &
                 -7* (iyear+(imonth+9)/12)       /4               &
                 -3*((iyear+(imonth-9)/ 7)/100+1)/4               &
               +275*         imonth              /9+iday+1721029
  !DaySecond
   DaT%Sec  = (iHour*60.d0+iMin)*60.d0+iSec*1.d0+imSec/1000.d0
  !Use DaySecond (& JulianDay) for differential timing to 1 millisec:
   return
End Subroutine Jdate20
!---------------------------------------------------------------------------7 9
