!read penmsh input files
!Generate inputfile for pentran
! main program
! complied by CVF 6.6 or ifc, pgf90 in windows
! complied by ifc 8.0 pgf90 v6 or higher on linux

!********************************************

!Main Program starting here
!******************************************************************
program inpred
use paraset1
use paraset2
use paraset4
use files
use ErrControl

ver_num="version 2.73b (June 2013)"
!*****************************************
!control varibles in this section

IsProcessFlux=0
IsProcessCurr=0
IsDotOut=0
IsOffTitan=0
!IsOffPentran=0

if( IsProcessFlux .eq. 1) then
 IsRefFlux=0
 IsCalPinPower=0
 IsPrbDotFlx(1)=0
 IsNomalizeFlux=0
endif

if(IsRefFlux .eq. 1) then
 IsPrbDotFlx(2)=0
endif

!IsSkipGeometry: 0, output geometry; 1, skip geometry output
IsSkipGeometry=0
!IsProcessFlux: 0, no; 1, read heart atn and act input files instead of penmsh input

! ****************************************
! define input files directory

input_dir='.'
flux_dir='.'
curr_dir='.'
ref_flux_dir='.'

!unit READLOG, logfile
outputfile(1)%fullname='read.log'

!*****************************************
!read input files

call GetArgument
call DisplayShortHeader

stime=SECNDS(0.0)

open(unit=READLOG, file=outputfile(1)%fullname)
write (READLOG,"('penmshxp ', A)" ) ver_num
write(form,"('read.log generated at : ')")
call TimeStamp(READLOG)

write(*,*)  'processing inputfiles.....'
if(IsProcessHrt .eq. 1) then
!in supplement01
 call ReadHrtInput

! stop 'give me a break, dude'
endif

if(IsProcessHrt .eq. 0) then
call ReadInput
 write(*,*)  'inputfiles done, logfile: ',trim(outputfile(1)%fullname)
 write(READLOG,*) 'PENMSH input files(.inp) completed'
 write(READLOG,*) '*************************************************'

 write(READLOG,*) 'Begin processing log'
 write(READLOG,*)

!***************************************
!processing input files in subcode2.f90
call Meshing
 write(*,*) 'Assigning Mat Num to fine meshes.....'
 call AssignMatNum
 write(*,*) 'Mat Num Assignment Done'
 write(READLOG,*) 'Assigning Mat Num to fine meshes.....done'
 if(s_format .ne. 2) then
   write(*,*) 'Assigning Src Mag. to fine meshes.....'
   call AssignSrcMag
   write(*,*) 'Src Mag Assignment Done'
   write(READLOG,*) 'Assigning Src Mag to fine meshes.....done'
 endif

endif
!***************************************
!processing flux if IsProcessFlux =1
if(IsProcessFlux .eq. 1) then
 write(*,*) 'allocating memory for flux...'
 call AllocateFluxArray
 write(*,*)  'allocation done'
 write(READLOG,*) 'allocating memory for flux...done'

 write(*,*) 'reading flux info...'
 call ReadInFluxset(1)
! call ReadInFlux
 write(*,*) 'flux done'
 write(READLOG,*) 'reading flux info...done'

 !avg. flux for each material
 if(IsFluxGM .ne. 0) call AvgMatFlux

 !flux output file
 if(IsFluxout .ne. 0 .and. num_group .eq. num_group_eff) &
  call MxOutputFlx
 !fov info for array project model
 if(IsFov .eq. 1) call WriteFov
 if(IsRefFlux .eq. 1) then
   write(*,*) 'reading Ref flux...'
   call ReadInFluxset(2)
   write(*,*) 'Ref flux done'
   write(READLOG,*) 'reading Ref. flux info...done'
 endif

endif

!***************************************
!processing current if IsProcessCurr =1
if(IsProcessCurr .eq. 1) then
 write(*,*) 'allocating memory for flux and current...'
 call AllocateFluxArray
 write(*,*)  'allocation done'
 write(READLOG,*) 'allocating memory for flux and current...done'

 write(*,*) 'reading flux&current info...'
 call ReadInCurrset(1)
! call ReadInFlux
 write(*,*) 'flux&current done'
 write(READLOG,*) 'reading flux&current info...done'
endif

!***************************************

!output geometry
if( IsSkipGeometry .eq. 0) then
  write(*,*) 'Dumping binary tecplot file...(-offplt option to turn off .plt file)'
  call MxOutputBinGemetry
  write(*,*) 'Geometry Done  ', trim(outputfile(3)%fullname)
  write(READLOG,*) 'Dumping binary tecplot file...done ', trim(outputfile(3)%fullname)
endif

!**************************************
! Write Dot f90 pentran input file
if (IsZSingle .eq. 0) then
if (IsProcessFlux .eq. 0) then

  if(IsOffPentran*IsOffTitan .eq. 0) then
    !assign transport parameters in subcode3.f90
     write(*,*) 'Assigning Transport Calculation Varibles...'
    call AssignPentranVaribles
    write(*,*) 'Assignment Done'
    write(READLOG,*) 'Assigning Pentran Varibles...done'

  endif

  if(IsOffPentran .eq. 0) then
    !Generate the pentran input file in subcode5.f90
    write(*,*) 'Generating pentran input deck...(-offf90 option to turn it off )'
    write(*,*) '(-nofido option to turn off fido in source section if it takes too long)'
    call WriteDotF90
    write(*,*) 'PENTRAN input file done ', trim(outputfile(4)%fullname)
    write(READLOG,*) 'Generating .pen pentran inputfile...done ',trim(outputfile(4)%fullname)
  else
    write(*,*) 'IsOffPentran=1: skipped generating PENTRAN input file'
    write(READLOG,*) 'skipped generating PENTRAN input file'
  endif

  if(IsOffTitan .eq. 0) then
    !Generate the TITAN input file in subcode5.f90
    write(*,*) 'Generating TITAN input deck...(-offtan option to turn off TITAN input )'
    call WriteTitanInp
    write(*,*) 'TITAN input file done: ', trim(outputfile(6)%fullname)
    write(READLOG,*) 'Generating TITAN input file...done ',trim(outputfile(6)%fullname)
  else
    write(*,*) 'IsOffTitan=1: skipped generating TITAN input file'
    write(READLOG,*) 'skipped generating TITAN input file'
  endif

elseif(IsProcessFlux .eq. 1) then
  write(*,*) 'skipped generating PENTRAN/TITAN input files'
  write(READLOG,*) 'skipped generating PENTRAN/TITAN input files'
endif
endif !IsZSingle plot

!output material .out file
if (IsDotOut .eq. 1) then
  write(*,*) 'generating .out file (penmsh outputfile) ...'
  call WriteDotOut
  write(*,*) '.out files done'
  write(READLOG,*) 'generating .out file (penmsh outputfile) ...done'
endif

if(IsPlotZlev .eq. 1) then
  write(*,*) 'Plotting z levels...(use -offpng option to turn off plotting)'
  call OutputPlot
  write(*,*) 'z lev graphic file(s)  done '
  write(READLOG,*) 'Plotting z level...  z lev graphic file(s) done '

endif

write(*,*)
call SubWarning
write(*,*)
write(*,"('All done, logfile : ',A )") trim(outputfile(1)%fullname)
write(READLOG,"('All done, logfile : ',A )") trim(outputfile(1)%fullname)
write(READLOG,*) '***************************************************'
stime=SECNDS(stime)
write(*,"('Total Runtime : ', f7.1, 'sec' )") stime
write(READLOG,"('Total Runtime : ', f7.1, 'sec' )") stime

write(READLOG,*)
write(form,"('run finished at : ')")
call TimeStamp(READLOG)
close(READLOG)

end

!***********************************************************************
!end of main flow

! command -i input_dir -f flux_dir_rel
! flux_dir=./input_dir/flux_dir_rel
subroutine GetArgument
use paraset1
use files
use funs
use ErrControl

parameter(ndash=2)
character*80,dimension(:), allocatable :: buf
integer count
integer i
integer :: IsSuccess=0
!dash_index :  flag     flag_argument
!      1        -i      input dir
!      2        -f      [rel flux dir]
!      3        -l      [logfile name]
integer,dimension(:),allocatable :: id_arg

real rarg
!count = iargc( )
count = IARGC( )
!buf(0): command itself
allocate(id_arg(0:count), buf(0:count))
id_arg=0
buf=''
call  GETARG(0,buf(0))

do i = 1, count

  call GETARG(i, buf(i))

  Flags: select case (trim(buf(i)))

    case ('-i','i','--input')
      id_arg(i)=1
    case ('-f', '-f1','-fl','--flux')
      id_arg(i)=2
      IsProcessFlux=1
    case ('-l','l','--log')
      id_arg(i)=3
    case ('-dot')
      id_arg(i)=4
    case('-w')
      id_arg(i)=5
    case('-nd','--ndmesh')
      id_arg(i)=6
	case ('-f2')
      id_arg(i)=7
      IsRefFlux=1
    case ('-f3', '-fjn')
     IsProcessCurr=1
     id_arg(i)=8
	case ('-n')
     IsNomalizeFlux=1
     id_arg(i)=9

	 case ('-ff')   !generate prbname.flx.out : flux summary
      IsFluxout=1
      IsProcessFlux=1
      id_arg(i)=11
    case ('-fa')   !generate prbname.adj.out : adjoint flux summary
      IsFluxout=2
      IsProcessFlux=1
      id_arg(i)=11

	case ('-chipath')   !generate prbname.adj.out : adjoint flux summary
      IsChiPath=1
      id_arg(i)=12

    case ('-offf90','-offpen')
      IsOffPentran=1
    case ('-offtan', '-offtitan')
      IsOffTitan=1
    case ('-prbflx','-prb')
      IsPrbDotFlx=1
    case('-nofido')
      IsFidoSrc=0

	case ('-offplt')
      IsSkipGeometry=1
    case ('-offpng')
      IsPlotZlev=0

	case ('-ming')
      flx_min_flag=2
    case ('-maxg')
      flx_max_flag=2

	case ('-min')
      id_arg(i)=28
      flx_min_flag=1
    case ('-max')
      id_arg(i)=29
      flx_max_flag=1

    case ('-plot', '-plotzyx')
      id_arg(i)=30
    case ('-xplot', '-plotx', '-xp')
      id_arg(i)=31
      IsXPlot=1
    case ('-yplot','-ploty', '-yp')
      id_arg(i)=32
      IsYPlot=1
    case ('-zplot','-plotz', '-zp')
      id_arg(i)=33
      IsZPlot=1
    case ('-pngsize','-png', '-size')
      id_arg(i)=34
    case('-plotmsf','-msf')
      id_arg(i)=35
    case ('-z')
      IsZSingle=1
      id_arg(i)=36
    case ('-color')
      id_arg(i)=37
    case ('-plot3d', '-3d','-3D')
      id_arg(i)=38
      Is3DPlot=1
    case ('-3dpos')
      id_arg(i)=39
      Is3DPlot=1
    case ('-nofm')
      IsDrawFm=0
    case('-log')
      IsLogPlot=1
    case ('-titan','-t')
      IsOffTitan=0
    case ('-gm', '-fgm')
      IsFluxGM=1
      IsProcessFlux=1
    case ('-agm')
      IsFluxGM=2
      IsProcessFlux=1
    case ('-fov')
      IsFov=1
 !for binary heart-like input
    case ('-hrt')
      IsProcessHrt=1
      id_arg(i)=40

    case ('-V', '-v')
      call ShowVersion

    case ('help','-h','-H','h','H','--help')
      call DisplayHelp
    case ('colormap','-helpcm','helpcm')
      call DisplayColorMap

    case default

      Flag_arg: select case (id_arg(i-1))

      case(1)  !-i
!         input_dir='./'//trim(buf(i))
         input_dir=trim(buf(i))
      case(2)  !-f1
         flux_dir=trim(buf(i))
      case(3)  !logfile name
         outputfile(1)%fullname=trim(buf(i))
      case(4)
         IsSuccess=CharNum(buf(i), DotBit)
         if(IsSuccess .eq. 0 ) then
            DotBit=-5
            write(*,"('GetArg: -dot flag argument fail','ignored -dot ', A,' using default Dotbit=-5 ')") &
              trim(buf(i))
         endif
      case(5)
         IsSuccess=CharNum(buf(i),max_warning)
         if(IsSuccess .eq. 0 ) then
            max_warning=5
            write(*,"('GetArg: -w flag argument fail, ignored -w ', A, ' using default max_warning=5')")&
              trim(buf(i))
         endif
      case(6)
         IsSuccess=CharNum(buf(i),ndmeth_global)
         if(IsSuccess .eq. 0 ) then
            ndmeth_global= 2
           write(*,"('GetArg: -nd flag argument fail, ignored -nd ',A, ' using default ndmeth_global=2 ')")&
              trim(buf(i))
        endif
      case(7) !-f2
         ref_flux_dir=trim(buf(i))
      case(8) !-f3
         curr_dir=trim(buf(i))
      case(9)
	    IsSuccess=CharFloat(buf(i),n_factor)
        if(IsSuccess .eq. 0 ) then
		  nfile%fullname=trim(buf(i))
		  n_factor=-1.0
		elseif (n_factor .lt. 0) then
          write(*,"('GetArg: -n flag argument fail, normalization factor has to be positive ',A, ' using default n_factor=max global flux ')")&
            trim(buf(i))
		  n_factor=0.0
		endif

	  case(11)
        IsSuccess=CharNum(buf(i),num_flx_out)
        if(IsSuccess .eq. 0 ) then
          num_flx_out=1
          write(*,"('GetArg: multi volumn flux file off ')")
        endif

	  case(12)
	    ChiPath_dir=trim(buf(i))
	  case(28)
	    IsSuccess=CharFloat(buf(i),flx_min)
		flx_min_flag=3
        if(IsSuccess .eq. 0 ) then
         write(*,"('GetArg: -min flag argument not present ',A, 'using default = min global flux')") &
           trim(buf(i))
		 flx_min_flag=1
!		elseif (flx_min .lt. 0) then
!          write(*,"('GetArg: -min flag argument fail, mininum flux has to be positive, ignored ',A)")&
!            trim(buf(i))
!		  flx_min=0.0
!		  flx_min_flag=1
		endif
      case(29)
	    IsSuccess=CharFloat(buf(i),flx_max)
		flx_max_flag=3
        if(IsSuccess .eq. 0 ) then
         write(*,"('GetArg: -max flag argument not present ',A,' using default = max global flux ')") &
           trim(buf(i))
		   flx_max_flag=1
!		elseif (flx_max .le. 0) then
!          write(*,"('GetArg: -max flag argument fail, maxium flux has to be positive, ignored ',A, ' using default n_factor=max global flux ')")&
!            trim(buf(i))
!		  flx_max=0.0
!		  flx_max_flag=1
		endif
	  case(30)
        IsSuccess=CharNum(buf(i),MidPlot)
        if(IsSuccess .eq. 0 ) then
          MidPlot=1
          write(*,"('GetArg: -plot flag argument fail, ignored ',A, ' using default plotting z levels ')")&
          trim(buf(i))
        endif
      case (31)
        call Str2Val(buf(i), rarg, IsSuccess)
        if(IsSuccess .eq. 0) then
          num_xloc=num_xloc+1
          id_arg(i)=31
          if(num_xloc .gt. maxplot) then
            write(*,"('GetArg: -xplot flag argument, too many positions, after ', I0, ' ignored')")   maxplot
          else
            xyzloc(num_xloc,1)=rarg
          endif
       else
          write(*,"('GetArg: -xplot flag argument, unregnized number: ',A, ' ignored')")   trim(buf(i))
       endif
     case (32)
       call Str2Val(buf(i), rarg, IsSuccess)
       if(IsSuccess .eq. 0) then
         num_yloc=num_yloc+1
         id_arg(i)=32
         if(num_yloc .gt. maxplot) then
           write(*,"('GetArg: -yplot flag argument, too many positions, after ', I0, ' ignored')")   maxplot
         else
           xyzloc(num_yloc,2)=rarg
         endif
       else
         write(*,"('GetArg: -yplot flag argument, unregnized number: ',A, ' ignored')")   trim(buf(i))
       endif
     case (33)
       call Str2Val(buf(i), rarg, IsSuccess)
       if(IsSuccess .eq. 0) then
         num_zloc=num_zloc+1
         id_arg(i)=33
         if(num_zloc .gt. maxplot) then
           write(*,"('GetArg: -zplot flag argument, too many positions, after ', I0, ' ignored')")   maxplot
         else
           xyzloc(num_zloc,3)=rarg
         endif
       else
          write(*,"('GetArg: -zplot flag argument, unregnized number: ',A, ' ignored')")   trim(buf(i))
       endif
    case(34)
      call Str2Val(buf(i), rarg, IsSuccess)
      if(IsSuccess .eq. 0) then
         if(rarg .gt. 9) then
            write(*,"('GetArg: -plotsize factor too big, using 9 ')")

         elseif(rarg .lt. 0.1) then
           write(*,"('GetArg: -plotsize factor too small , using 0.1 ')")
         else
           PngSizeFactor=rarg
         endif
      else
        write(*,"('GetArg: -pngsize flag argument, unregnized number: ',A, ' ignored')")   trim(buf(i))
      endif
    case (35)
     IsSuccess=CharNum(buf(i),MatSrcFlx)
     if(IsSuccess .eq. 0 ) then
       MatSrcFlx=1
       write(*,"('GetArg: -plot flag argument fail, ignored ',A, ' using default plotting z levels ')")&
       trim(buf(i))
     endif
     if (abs(MatSrcFlx) .ge. 4) IsProcessFlux=1
   case (36)
       IsSuccess=CharNum(buf(i),z_start)
       z_end=z_start
       if(IsSuccess .eq. 0 ) then
          z_start=1
          z_end=1
          write(*,"('GetArg: -z option without an argument, plotting the frist z level only ')")
       endif
   case(37)
     IsSuccess=CharNum(buf(i),ColorMap)
     if(IsSuccess .eq. 0 ) then
        ColorMap=3
        write(*,"('GetArg: invalid colormap id, use <penmsh colormap> for more help, set to default')")
     endif
   case (38)
      IsSuccess=CharFloat(buf(i),view_ang)
	  if(Issuccess .eq. 0) then
	  view_ang=0
	  write(*,"('GetArg: -3d option without an argument, 3d plotting from 0 degree angle (default) ')")
      endif
   case (39)
    call Str2Val(buf(i), rarg, IsSuccess)
	if(IsSuccess .eq. 0) then
      num_vpos=num_vpos+1
      id_arg(i)=39
      if(num_vpos .gt. 3) then
        write(*,"('GetArg: -3dpos flag argument, too many values, first three entries are used')")
      else
       view_pos(num_vpos)=rarg
      endif
    else
      write(*,"('GetArg: -3dpos without argument, using view point (3,3,3) default')")
    endif

   case (40)
       hrt_filename=trim(buf(i))
   case default
         write(*,*) "Warning: ", trim(buf(i)), ": unknown argument ignored"
         write(*,*) "Use inpred -h to display help"
         stop
   end select Flag_arg

  end select  Flags

enddo

if (IsProcessFlux .eq. 0 .and. IsRefFlux .eq. 1) then
  write(*,*) "Commandline option conflict: ",  "-f2  option used only if -f1 presents"
  write(*,*) "Use inpred -h to display help"
  stop
endif

if (IsProcessFlux .eq. 0 .and. IsRefFlux .eq. 1) then
  write(*,*) "Commandline option conflict: ",  "-f2  option used only if -f1 presents"
  write(*,*) "Use inpred -h to display help"
  stop
endif

if (IsZSingle .eq. 1) then

if(IsProcessHrt .eq. 1) then
  write(*,*) "Commandline option conflict: ",  "-z  and -hrt'"
  write(*,*) "Use inpred -h to display help"
  stop
endif

endif

If(IsProcessCurr .eq. 1) then
 if(IsProcessFlux .ne. 0 .or. IsRefFlux .ne. 0) then
  write(*,*) "Commandline option conflict: ",  "-f3 (fjn) option is used, -f1 or -f2 will be ignored "
  write(*,*) "Use inpred -h to display help"
  IsProcessFlux=0
  IsRefFlux=0
 endif
endif

end subroutine

!Display the help screen and terminate the program
subroutine DisplayHelp
use paraset1, only: ver_num

write(*,"(A)") "NAME"
write(*,"(A)")
write(*,"(A)") "PENMSHXP - generate PENTRAN inputfile"
write(*,'("           (inpred)  " ,A)')  trim(ver_num)
write(*,"(A)") "   input files : penmsh.inp, prb_name[z_level#].inp"
write(*,"(A)") "      optional : prb_name[grp#].flx or prb[grp#].flx : "
write(*,"(A)") "                    pentran serial/parallel  flux files "
write(*,"(A)") "                 prb_name.spc, prb_name.chi   : "
write(*,"(A)") "                    src spectrum and fission spectrum file "
write(*,"(A)") "                 prb_name.mba (mat balance file) "
write(*,"(A)") "                 prb_name.src (src density grid file) "
write(*,"(A)")
write(*,"(A)") "  output files : prb_name_out.pen (PENTRAN input); "
write(*,"(A)") "                 prb_name_mix.plt (tecplot data file) "
write(*,"(A)") "                 prb_name.mcr (tecplot macro file) "
write(*,"(A)") "                 prb_name_out.mba (Material vol balance ) "
write(*,"(A)") "                 read.log         (read log file        ) "
write(*,"(A)")
write(*,"(A)") "SYNOPSIS"
write(*,"(A)")
write(*,"(A)") "    penmshxp [flag1] [flag1 argument] [flag2] [flag2 argument] ..."
write(*,"(A)")
write(*,"(A)") "DESCRIPTION"
write(*,"(A)")
write(*,"(A)")  'Flag          Flag_argument'
write(*,"(A)")
write(*,"(A)")  '-i,--input    input_dir'
write(*,"(A)")  '       input file directory, default current dir '
write(*,"(A)")
write(*,"(A)")  '-f1,-f,--flux     flux_dir'
write(*,"(A)")  '       flux file dir, default: current dir '
write(*,"(A)")  '       ./[flux_dir]'
write(*,"(A)")  '-f2           ref_flux_dir'
write(*,"(A)")  '       reference flux file dir, default: current dir '
write(*,"(A)")  '       ./[flux_dir]'
write(*,"(A)")  '-f3, -fjn         curr_dir'
write(*,"(A)")  '       flux current file dir, default: current dir '
write(*,"(A)")  '       ./[flux_dir]'
write(*,"(A)")  ' Note: .fjn files are geneated by PENDATA (Option 9) with 5 rows only:  '
write(*,"(A)")  '       Phi   J    Jx    Jy     Jz (Select Varibles No.10-14 in PENDATA) '
write(*,"(A)")  '-n     <factor or a filename  > '
write(*,"(A)")  '   To nomalized fluxes flux=flux/factor, if no argument given, normalized to global max flux'
write(*,"(A)")  '   filename: a file containing a factor array for each group (fido supported) '
write(*,"(A)")  '             size=num_group, in increasing order (fwd) or decreasing order (adj) '
write(*,"(A)")  '   e.g. -n 100.0  all fluxes are divided by 100,'
write(*,"(A)")
write(*,"(A)")  '-l,--log      logfile name'
write(*,"(A)")  '       log file name'
write(*,"(A)")  '-dot   DotBit   '
write(*,"(A)")  '       FIDO control on spacef value format'
write(*,"(A)")  '       DotBit=3, i.e 45.123, 333.123   '
write(*,"(A)")  '       DotBit=-5 i.e  1.12345E-01 ,default        '
write(*,"(A)")  '-nd --ndmeth  ndmeth   '
write(*,"(A)")  '       GLOBLE ndmeth varible: differencing scheme '
write(*,"(A)")  '       ndmeth=2 : default, adaptive DTW, upgradable to EDW   '
write(*,"(A)")  '       ndmeth=-2  fixed DTW       '
write(*,"(A)")  '-offplt                         '
write(*,"(A)")  '     Skip tecplot Geometry file  '
write(*,"(A)")  '-offpng                         '
write(*,"(A)")  '     Skip z level plots   '
write(*,"(A)")  '-plotmsf -msf MatSrcFlx                        '
write(*,"(A)")  '    MatSrcFlx=1, 2 or 4: plotting mat, src, or flx  '
write(*,"(A)")  '    default=1 , e.g. -plotmsf 5 will plot mat and flx  '
write(*,"(A)")  '    e.g -plotmsf -4 will plot flux, and write flux to a .cat file  '
write(*,"(A)")  '-plotzyx MidPlot                        '
write(*,"(A)")  '    MidPlot=1,2 or 4: plotting z, y, or x mid-levels  '
write(*,"(A)")  '    default=1 , e.g. -plotzyx 3 will plot z and y  '
write(*,"(A)")  '-plotx  pos1 pos2...                         '
write(*,"(A)")  '    x-y plot only at z= <pos1 pos2 ...>     '
write(*,"(A)")  '-ploty  pos1 pos2 ...                         '
write(*,"(A)")  '    x-z plot only at y= <pos1 pos2 ...>     '
write(*,"(A)")  '-plotz  pos1 pos2...                         '
write(*,"(A)")  '    y-z plot only at x= <pos1 pos2 ...>     '
write(*,"(A)")  '-size  factor                         '
write(*,"(A)")  '    png file resolution multiplication factor     '
write(*,"(A)")  '-3d  view_angle                         '
write(*,"(A)")  '    generate 3d plots, and set view angle in degree     '
write(*,"(A)")  '    default angle is 0 degree and top view, negative entry for bottom view  '
write(*,"(A)")  '-3dpos  x  y  z                         '
write(*,"(A)")  '    generate 3d plots, and set viewpoint position at (x, y, z) '
write(*,"(A)")  '    default position is (4 4 4),  plotting area coordinates  '

write(*,"(A)")  '-max  <flux max>                         '
write(*,"(A)")  '    flux plot scale : maxium, default max among all groups   '

write(*,"(A)")  '-min  <flux min>                         '
write(*,"(A)")  '    flux plot scale : minmum, warning lower than the minmum points will black out '

write(*,"(A)")  '-maxg                          '
write(*,"(A)")  '    flux plot scale : max for one group   '

write(*,"(A)")  '-ming                           '
write(*,"(A)")  '    flux plot scale : minmum for one group '

write(*,"(A)")  '-nofm                         '
write(*,"(A)")  '    turn off Fine Mesh lines in all coarse meshes   '
write(*,"(A)")  '    Note: FM lines automatically off if FM size < 8 pixels in a CM '
write(*,"(A)")  '-offf90                         '
write(*,"(A)")  '    Skip generate pentran input file  '
write(*,"(A)")  '-z  z_level_num                       '
write(*,"(A)")  '    plot one z-level only  '
write(*,"(A)")  '-color ColorMap number'
write(*,"(A)")  '   type penmshxp -helpcm for more help'

write(*,"(A)")  '-fgm                         '
write(*,"(A)")  '    generate file flux.fgm : avg. flux per material zone for each group'
write(*,"(A)")  '    Note: require load the flux using -f option '
write(*,"(A)")  '-agm                         '
write(*,"(A)")  '    generate file flux.agm: avg. adjoint flux per mat zone for each group'
write(*,"(A)")  '    Note: groups are flipped  '
write(*,"(A)")  '-ff                         '
write(*,"(A)")  '    generate file prbname.flx.out : flux distribution file'
write(*,"(A)")  '    Note: require load the all group flux using -f option'
write(*,"(A)")  '-fa                         '
write(*,"(A)")  '    generate file prbname.adj.out : adjoint flux distribution file '
write(*,"(A)")  '    Note: groups are flipped  '

write(*,"(A)")  '-titan                         '
write(*,"(A)")  '     generate titan input file  '
!write(*,"(A)")  '-offtitan                         '
!write(*,"(A)")  '     Skip generate titan input file  '
write(*,"(A)")  '-nofido                         '
write(*,"(A)")  '     Turn off FIDO generation in source section of PENTRAN input deck '
write(*,"(A)")  '-w   number                     '
write(*,"(A)")  '     Max warnings                '
write(*,"(A)")  '-hrt   phantom input file name    '
write(*,"(A)")  '    handle binary phantom data files                 '
write(*,"(A)")  '-V       '
write(*,"(A)")  '    show version number                 '
write(*,"(A)")

stop  'Report bugs to CE YI<yice@gatech.edu>'

end subroutine

subroutine DisplayColorMap

write(*,"(A)") 'Use -color <colormap id> option to set a color map'
write(*,"(A)") 'colormap id: '
write(*,"(A)") ' 0 : dislin id=SMALL, defines a small colour table with the 8 colours'
write(*,"(A)") '     1 = BLACK, 2 = RED, 3 = GREEN, 4 = BLUE, 5 = YELLOW, 6 = ORANGE, 7 = CYAN and 8 = MAGENTA'
write(*,"(A)") ' 1 : dislin id= VGA, defines the 16 standard colours of a VGA graphics card'
write(*,"(A)") ' 2 : dislin id=RAIN, defines 256 colours arranged in a rainbow where 0 means black and 255 means white'
write(*,"(A)") ' 3 : dislin id=SPEC, defines 256 colours arranged in a rainbow where 0 means black and 255 means white'
write(*,"(A)") '     This colour table uses more violet colours than RAIN'
write(*,"(A)") ' 4 : dislin id=GREY, defines 256 grey scale colours where 0 means black and 255 is white'
write(*,"(A)") ' 5 : dislin id=RRAIN, is the reverse colour table of RAIN '
write(*,"(A)") ' 6 : dislin id=RSPEC, is the reverse colour table of SPEC'
write(*,"(A)") ' 7 : dislin id=RGREY, is the reverse colour table of GREY'
write(*,"(A)") ' 8 : dislin id=TEMP, defines a temperature colour table'
write(*,"(A)") ' The default colour table is 3 (SPEC) '
write(*,"(A)")

stop  'Report bugs to CE YI<yice@gatech.edu>'

return
end subroutine
!Display the help screen and terminate the program
subroutine DisplayHeader
use paraset1, only: ver_num

write(*,"(A)")
write(*,"(A)")
write(*,"(A)") "                              PENMSHXP                             "
write(*,"(A,A)") "                PENMSH Express ", trim(ver_num)
write(*,"(A)") "           A Mesh Generator to Build PENTRAN Input Deck            "
write(*,"(A)") "     with compatibility to PENMSH developed by Dr. A. Haghighat    "
write(*,"(A)") "            (DISLIN graphic library by Helmut Michels)             "
write(*,"(A)")
write(*,"(A)")
write(*,"(A)") "                               CE YI                              "
write(*,"(A)") "                        email: yice@gatech.edu                    "
write(*,"(A)")
write(*,"(A)")
write(*,"(A)") "                    Use 'penmshxp -h' for more info               "
write(*,"(A)")
write(*,"(A)")
write(*,"(A)")
write(*,"(A)")
write(*,"(A)")  'Begin processing...'

end subroutine

!Display the help screen and terminate the program
subroutine DisplayShortHeader
use paraset1, only: ver_num

write(*,"('PENMSH ', A)")  trim(ver_num)
write(*,"('use -h option for help, and check the read.log file')")
write(*,"(A)")
write(*,"(A)")  'Begin processing...'

end subroutine

!Display Version and copyright info
Subroutine ShowVersion
use paraset1, only : ver_num

write(*, "(' PENMSHXP ', A)")  ver_num
write(*, "(' Copyright 2011 Ce Yi at GaTech')" )
if( index(ver_num, 'b ') .ne. 0)   &
write(*, "(' DISLIN Graphic Library by Helmut Michels at Max Planck Institute' )")
stop '      '
end subroutine

!put a tiem stamp on unit_num
!if unit_num=0, stamp on the screen
!a message delivered too at the same time
!the message is stored in 'form' @ module 'files'
subroutine TimeStamp(unit_num)
use files

integer unit_num

character(LEN =8) str_date
character(LEN =12) str_time

call date_and_time (str_date, str_time)

if(unit_num .eq. 0) then

write(*, "(A, 4x,A2,':',A2,':',A2,2x,'on ',A2,'/',A2,'/', A4)" )  trim(form), &
   str_time(1:2),str_time(3:4),str_time(5:6),str_date(5:6),str_date(7:8), str_date(1:4)
write(*,*)

else
write(unit_num, "(A, 4x,A2,':',A2,':',A2,2x,'on ',A2,'/',A2,'/', A4)" )  trim(form), &
   str_time(1:2),str_time(3:4),str_time(5:6),str_date(5:6),str_date(7:8), str_date(1:4)
write(unit_num,*)
endif

end subroutine

!convert string to value
subroutine Str2Val(str, val, ierr)

use ErrControl

character(*) str
real val
integer ierr

read(str,*, iostat=ierr) val
!if(ierr .ne. 0) then
!write(err_message,"('error in converting string to val ')")

!call TrapInputError(1)
!endif

return
end subroutine


text:
*
*             T R I P L E
*             ***********
*
*
*       Three-body regularization program.
*       ----------------------------------
*
*       Method of Aarseth & Zare, Celestial Mechanics 10, 185.
*       ......................................................
*
*       Developed by Sverre Aarseth, IOA, Cambridge.
*       ............................................
*
      FUNCTION AarsethZare(TIME,BODY,POS,VEL,TCRIT,Etot,NREG)
c      FUNCTION AarsethZare(TIME,M,X,XDOT)
c      PROGRAM TRIPLE
*
      IMPLICIT  REAL*8  (A-H,M,O-Z)
      REAL*8 TIME, body(3), pos(3,3), vel(3,3), Etot, TCRIT
      INTEGER NREG
      COMMON/AZREG/  Q(8),P(8),R,R1,R2,ENERGY,M(3),X(3,3),XDOT(3,3),
     &               RCOLL,ERROR,C11,C12,C19,C20,C24,C25,NSTEPS,NAME(3)
      COMMON/CLOSE/  RIJ(3,3),ICALL
      COMMON/BSSAVE/  EP(4),TFAC,ITFAC,JC,NHALF2
      COMMON/NCALL/  NFN
      REAL*8  Y(17),ORDER(16),RK(3)
      INTEGER AarsethZare

      write (*,*) "Input:", TIME,BODY,POS,VEL,TCRIT,Etot,NREG
*
      AarsethZare = 1
      Etot = 0
*
*       COMMON variables
*       ****************
*
*       ------------------------------------------------------------------
*       C11     Inverse mass factor for DERQP (also C12,C19,C20,C24,C25).
*       ENERGY  Twice the initial total energy.
*       ERROR   Relative energy error at the end (ignore with T' = 1/U).
*       ICALL   Indicator for pericentre check (first function call only).
*       M       Particle mass.
*       NAME    Particle identity (initialized to 1,2,3).
*       NSTEPS  Number of DIFSY calls.
*       P       Regularized momenta.
*       Q       Regularized coordinates.
*       R       Distance between M(1) and M(2) (not regularized).
*       R1      Distance between M(1) and M(3).
*       R2      Distance between M(2) and M(3).
*       RCOLL   Minimum two-body separation (osculating pericentre).
*       RIJ     Minimum pairwise separations (osculating pericentre).
*       X       Particle coordinates (X(3,I) is Z-component).
*       XDOT    Velocity components (XDOT(3,I) is Z-component).
*       ------------------------------------------------------------------
*
*
      NEXP = 0
*       Read tolerance and termination time.
      TOL0 = 1.e-12
cccc  TCRIT = 100
c      READ (5,*)  TOL0, TCRIT
*
*       Generate next initial condition (place most active body in M(3)).
c      DO 2 I = 1,3
c          M(I) = body(i)
c          DO 3 K = 1,3
c             X(K,I) = pos(k,i)
c             XDOT(K,I) = vel(k,i)
c    3     CONTINUE
c    2 CONTINUE

      CALL DATA(body, pos, vel)
*
*       Save the tolerance.
      EPS = TOL0
*
*       Initialize diagnostic variables & counters.
      R12MIN = 100.0
      RMIN = 100.0
      RCOLL = 100.0
      DO 10 J = 1,3
          DO 5 K = 1,3
              RIJ(J,K) = 100.0
    5     CONTINUE
   10 CONTINUE
      NSTEPS = 0
      NREG = 0
      ICALL = 0
      NFN = 0
      JC = -1
*
*       Initialize local time & regularized time.
      TIME = 0.0D0
      TAU = 0.0D0
      Y(17) = 0.0D0
*       Specify the number of first-order equations for the integrator.
      N = 17
*
*       Obtain initial energy and transform to regularized variables.
      CALL TRANSF(1)
*
*       Define gravitational radius and pericentre check distance.
      RGRAV = (M(1)*M(2) + M(1)*M(3) + M(2)*M(3))/(0.5D0*ABS(ENERGY))
      RSTAR = 0.5*RGRAV
*
*       Form the two smallest distances (assume sensible reference body).
      R1 = Q(1)**2 + Q(2)**2 + Q(3)**2 + Q(4)**2
      R2 = Q(5)**2 + Q(6)**2 + Q(7)**2 + Q(8)**2
*
*       Set termination distance to maximum of separation and 10*RGRAV.
      RMAX0 = 1.1*MAX(R1,R2,10.0*RGRAV)
*
*       Specify the crossing time (also meaningful if ENERGY > 0).
      TCR = (M(1) + M(2) + M(3))**2.5/ABS(ENERGY)**1.5
      RAP = 10.0*RMAX0
      TNEXT = 0.0
*
*       Define a nominal crossing time for nearly parabolic systems.
      RMAX = MAX(R1,R2)
      TSTAR = RMAX*SQRT(RMAX/(M(1) + M(2) + M(3)))
*
*       Determine the smallest two-body time-scale from parabolic orbit.
      IM = 1
      RM = R1
      IF (R2.LT.R1) THEN
          IM = 2
          RM = R2
      END IF
      VP2 = 2.0*(M(IM) + M(3))/RM
      TP = RM/SQRT(VP2)
      TSTAR = MIN(TP,TSTAR)
c      write (*,*) "TIMES=", TP, TSTAR
      TORB = TP
*     TCRIT = TORB*TCRIT
c     Reinstated: integrate for 100 orbits
c     Assure that TCRIT is properly defined if not given as input (Sept.2010).
      if (TCRIT.le.0) then
         TCRIT = TORB*100.0
      endif

*       Set TPR and initial step (ignore M1*M2/R if DT/DTAU = 1/U).
*     TPR = R1*R2/(M(1)*M(3)*R2 + M(2)*M(3)*R1)
      TPR = R1*R2/SQRT(R1 + R2)
      DTAU = MIN(TCR,TSTAR)*EPS**0.1/TPR
*
*       Try alternative expression for initial step (Seppo Mikkola).
*     rmxt=max(r1,r2)
*     dtau=.1*(EPS/1.0E-12)**0.1*rmxt*sqrt(rmxt/(m(1)+m(2)+m(3)) )/tpr

*       Initialize time constant & input array for the integrator.
      CONST = 0.0D0
      DO 20 K = 1,8
          CONST = CONST + Q(K)*P(K)
          Y(K) = Q(K)
          Y(K+8) = P(K)
   20 CONTINUE
*
*       Produce initial output.
      CALL TRANSF(2)
*
*       Specify reference values for each vector (Mikkola's method).
      DO 35 K = 1,4
          K1 = 4*(K - 1)
          SN1 = 0.25*(ABS(Y(K1+1)) + ABS(Y(K1+2)) + ABS(Y(K1+3)) +
     &                                              ABS(Y(K1+4)))
*       Initialize the reference vector.
          IF (NSTEPS.EQ.0) THEN
              ORDER(K1+1) = SN1
          END IF
*       Do not permit any small values because of relative tolerance.
          SN1 = 0.1*SN1 + 0.9*ORDER(K1+1)
          DO 34 L = 1,4
              ORDER(K1+L) = SN1
   34     CONTINUE
   35 CONTINUE
*
*       Advance the equations of motion by Bulirsch-Stoer integrator.
*     CALL DIFSY1(N,EPS,ORDER,DTAU,TAU,Y)
   30 CALL DIFSY1(N,EPS,DTAU,TAU,Y)
*
*       Copy regularized coordinates & momenta and obtain physical time.
c      SUMQP = 0.0D0
      DO 40 K = 1,8
          Q(K) = Y(K)
          P(K) = Y(K+8)
c          SUMQP = SUMQP + Q(K)*P(K)
*       Note that the momentum includes factor of 2 in AZ formulation.
   40 CONTINUE
*
*       Set explicit time (Baumgarte & Stiefel, 1974 & Aarseth, 1976).
c      TIME = (0.5D0*(SUMQP - CONST) - TAU)/ENERGY
       TIME = Y(17)
      write (*,*) "TIMES=", TP, TSTAR, TIME, X

*
*       Update relative distances (NB! not quite latest value).
      R1 = Q(1)**2 + Q(2)**2 + Q(3)**2 + Q(4)**2
      R2 = Q(5)**2 + Q(6)**2 + Q(7)**2 + Q(8)**2
*
*       Check minimum two-body separations and increase step counter.
      RMIN = MIN(RMIN,R)
      RM = MIN(R1,R2)
      R12MIN = MIN(R12MIN,RM)
      RMAX = MAX(R1,R2,R)
      NSTEPS = NSTEPS + 1
*
*       Check minimum two-body separations.
      RK(1) = R1
      RK(2) = R2
      RK(3) = R
*       Consider pairs 1-2, 1-3 & 2-3 with identified names.
      DO 44 K = 1,3
          DO 42 L = K+1,3
              I = NAME(K)
              J = NAME(L)
*       Use cyclic loop index (3,1,2) for distances R, R1 & R2.
              KK = K - 1
              IF (KK.EQ.0) KK = 3
              RIJ(I,J) = MIN(RIJ(I,J),RK(KK))
              RIJ(J,I) = MIN(RIJ(J,I),RK(KK))
   42     CONTINUE
   44 CONTINUE
*
*       Switch on search indicator inside RSTAR (reset in DERQP).
      IF (RM.LT.RSTAR) THEN
          ICALL = 1
      END IF
*
*       Use actual two-body force to decide on switching.
      F12 = (M(1) + M(2))/R**2
      F13 = (M(1) + M(3))/R1**2
      F23 = (M(2) + M(3))/R2**2
      IF (F12.LT.MAX(F13,F23)) GO TO 70
*
      IMIN = 1
*       Use a simple distance test to determine new reference body IMIN.
      IF (R2.LT.1.00001*R1) IMIN = 2
*
*       Transform to physical variables and rename the exchanged particles.
      CALL TRANSF(3)
*
      DO 50 K = 1,3
          TEMP1 = X(K,3)
          TEMP2 = XDOT(K,3)
          X(K,3) = X(K,IMIN)
          XDOT(K,3) = XDOT(K,IMIN)
          X(K,IMIN) = TEMP1
          XDOT(K,IMIN) = TEMP2
   50 CONTINUE
*
      TEMP1 = M(3)
      M(3) = M(IMIN)
      M(IMIN) = TEMP1
      NAME3 = NAME(3)
      NAME(3) = NAME(IMIN)
      NAME(IMIN) = NAME3
*
*       Transform back to regularized variables and initialize input array.
      CALL TRANSF(4)
      DO 60 K = 1,8
          Y(K) = Q(K)
          Y(K+8) = P(K)
   60 CONTINUE
*
*       Update regularization counter at the end of switching procedure.
      NREG = NREG + 1
*
*       Check termination criteria (TIME > TCRIT or RMAX > RMAX0).
   70 IF (TIME.GT.TCRIT) GO TO 90
*     IF (MAX(R1,R2).LT.RAP) GO TO 30
*     IF (TIME.LT.TNEXT) GO TO 30
      TNEXT = TNEXT + TORB
*
*       Obtain final output after transforming to physical variables.
      CALL TRANSF(2)
*
*       See whether the final configuration is well resolved.
      RATIO = R1/R2
*      IF (RATIO.GT.0.1.AND.RATIO.LT.10.0) GO TO 90
      IF (RATIO.GT.0.1.AND.RATIO.LT.10.0) GO TO 30
*
*       Set index of second binary component & escaper.
      IMIN = 1
      IF (R2.LT.R1) IMIN = 2
      IESC = 3 - IMIN
*
*       Evaluate orbital elements.
      RDOT = 0.0D0
      VREL2 = 0.0D0
      VESC = 0.0D0
      DO 75 K = 1,3
          RDOT = RDOT + (X(K,3) - X(K,IMIN))*(XDOT(K,3) - XDOT(K,IMIN))
          VREL2 = VREL2 + (XDOT(K,3) - XDOT(K,IMIN))**2
          VESC = VESC + (XDOT(K,3) - XDOT(K,IESC))**2
   75 CONTINUE
*
      RB = RM
      MB = M(3) + M(IMIN)
      SEMI = 2.0D0/RB - VREL2/MB
      VESC = SQRT(VESC)/SQRT(MB*ABS(SEMI))
*       Velocity of escaper w.r. to M(3) (scaled by binary velocity).
      SEMI = 1.0/SEMI
      E = SQRT((1.0 - RB/SEMI)**2 + RDOT**2/(SEMI*MB))
      EB = -0.5D0*M(3)*M(IMIN)/SEMI
      IF (TIME.LT.TORB) SEMI0 = SEMI
*
      IF (EB.LT.0.5*ENERGY) THEN
          EBE = EB/(0.5*ENERGY)
*       Binding energy scaled by total energy.
*
          WRITE (6,80)  TIME, NAME(IMIN), NAME(3), R1, R2,
     &                  SEMI, E, NREG, ERROR, Etot
   80     FORMAT  ('TIME= ', E12.4, ' BINARY ',2I3,'  R1 =',1P,E9.1,
     &                  '  R2 =',E9.1,
     &                  '  A =',E12.4,'  E =',0P,F7.3,'  NREG =',I3,
     &                  '  DE/E =',E9.1,'  DEtot/Etot =',E9.1)
c
c         copy the final parameters
          DO 85 L = 1,3
             I = NAME(L)
             BODY(I) = M(L)
             DO 84 K = 1,3
                POS(K,I) = X(K,L)
                VEL(K,I) = XDOT(K,L)
   84        CONTINUE
   85     CONTINUE

c         Last two lines added
          AarsethZare = 0
          return
          GO TO 90
      END IF
*     IF (TIME.LT.TCRIT.AND.ABS(SEMI-SEMI0).LT.0.5*SEMI0) GO TO 30
      IF (TIME.LT.TCRIT) GO TO 30
*
   90 NEXP = NEXP + 1
      CALL TRANSF(2)
*
c
c
*       Identify particles and copy  mass, coordinates and velocities.
      DO 95 L = 1,3
          I = NAME(L)
          BODY(I) = M(L)
          DO 94 K = 1,3
              POS(K,I) = X(K,L)
              VEL(K,I) = XDOT(K,L)
   94     CONTINUE
   95 CONTINUE
c
c
      RCOLL = MIN(RCOLL,R12MIN)
      Etot = Etot + ERROR
      WRITE (*,*)  "TORB+", TIME, TORB, RIJ(1,2), RIJ(1,3), RIJ(2,3)
      WRITE (6,100)  TIME/TORB, TIME, RMIN, R12MIN, ERROR, NSTEPS
  100 FORMAT  (/,'  T/TK =',F7.1, ' TIME=', F7.1, ' MIN(R) =',1PE8.1,
     &        '  MIN(R1,R2) =',E8.1,
     &                 '  DE/E =',E9.1,'  NSTEPS =',I8)
c      WRITE (6,105)  RIJ(1,2), RIJ(1,3), RIJ(2,3), RCOLL
c  105 FORMAT (' RIJ:  1-2 1-3 2-3 RCOLL ',1P4E10.2)
c      WRITE (6,110)  NSTEPS, NFN
c  110 FORMAT ('  NSTEPS =',I9,'  NFN =',I11)
*
*       Continue with the next experiment.
c      GO TO 1
      AarsethZare = 0
      return

*
      END
      include 'difsy1.f'
      include 'data.f'
      include 'peri.f'
      include 'block.f'
      include 'derqp.f'
      include 'transf.f'

ext: f

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     https://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.

module test_curve

  use, intrinsic :: iso_c_binding, only: c_bool, c_double, c_int
  use curve, only: &
       CurveData, evaluate_curve_barycentric, evaluate_multi, &
       specialize_curve, evaluate_hodograph, subdivide_nodes, newton_refine, &
       LOCATE_MISS, LOCATE_INVALID, locate_point, elevate_nodes, &
       get_curvature, reduce_pseudo_inverse, full_reduce, compute_length, &
       curves_equal, subdivide_curve
  use types, only: dp
  use unit_test_helpers, only: &
       MACHINE_EPS, print_status, get_random_nodes, get_id_mat
  implicit none
  private &
       test_evaluate_curve_barycentric, test_evaluate_multi, &
       test_specialize_curve, test_evaluate_hodograph, test_subdivide_nodes, &
       subdivide_points_check, test_newton_refine, test_locate_point, &
       test_elevate_nodes, test_get_curvature, test_reduce_pseudo_inverse, &
       pseudo_inverse_helper, test_full_reduce, test_compute_length, &
       test_curves_equal, test_subdivide_curve
  public curve_all_tests

contains

  subroutine curve_all_tests(success)
    logical(c_bool), intent(inout) :: success

    call test_evaluate_curve_barycentric(success)
    call test_evaluate_multi(success)
    call test_specialize_curve(success)
    call test_evaluate_hodograph(success)
    call test_subdivide_nodes(success)
    call test_newton_refine(success)
    call test_locate_point(success)
    call test_elevate_nodes(success)
    call test_get_curvature(success)
    call test_reduce_pseudo_inverse(success)
    call test_full_reduce(success)
    call test_compute_length(success)
    call test_curves_equal(success)
    call test_subdivide_curve(success)

  end subroutine curve_all_tests

  subroutine test_evaluate_curve_barycentric(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes(3, 4)
    real(c_double) :: lambda1(3), lambda2(3)
    real(c_double) :: evaluated(3, 3), expected(3, 3)
    integer :: case_id
    character(26) :: name

    case_id = 1
    name = "evaluate_curve_barycentric"

    ! CASE 1: Slightly complex example.
    nodes(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
    nodes(:, 2) = [0.5_dp, 3.0_dp, 1.0_dp]
    nodes(:, 3) = [1.5_dp, 4.0_dp, 1.0_dp]
    nodes(:, 4) = [2.0_dp, 8.0_dp, 1.0_dp]
    lambda1 = [0.25_dp, 0.5_dp, 0.75_dp]
    lambda2 = [0.25_dp, 0.125_dp, -0.75_dp]
    expected(:, 1) = [0.125_dp, 0.453125_dp, 0.109375_dp]
    expected(:, 2) = [0.0859375_dp, 0.390625_dp, 0.119140625_dp]
    expected(:, 3) = [0.421875_dp, -2.109375_dp, -0.421875_dp]
    call evaluate_curve_barycentric( &
         4, 3, nodes, 3, lambda1, lambda2, evaluated)

    case_success = all(evaluated == expected)
    call print_status(name, case_id, case_success, success)

  end subroutine test_evaluate_curve_barycentric

  subroutine test_evaluate_multi(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(3, 2)
    real(c_double) :: s_vals1(129)
    real(c_double) :: evaluated1(3, 129), expected1(3, 129)
    real(c_double) :: nodes2(2, 3)
    real(c_double) :: s_vals2(65)
    real(c_double) :: evaluated2(2, 65), expected2(2, 65)
    integer :: i
    integer :: case_id
    character(14) :: name

    case_id = 1
    name = "evaluate_multi"

    ! CASE 1: Linear curve.
    do i = 1, 129
       s_vals1(i) = (i - 1) / 128.0_dp
    end do
    nodes1(:, 1) = [1.0_dp, 1.0_dp, -7.0_dp]
    nodes1(:, 2) = [2.0_dp, -1.0_dp, -4.0_dp]
    ! B(s) = [s + 1, 1 - 2 s, 3 s - 7]
    expected1(1, :) = 1.0_dp + s_vals1
    expected1(2, :) = 1.0_dp - 2.0_dp * s_vals1
    expected1(3, :) = -7.0_dp + 3.0_dp * s_vals1
    call evaluate_multi( &
         2, 3, nodes1, 129, s_vals1, evaluated1)
    case_success = all(evaluated1 == expected1)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic curve.
    do i = 1, 65
       s_vals2(i) = (i - 1) / 64.0_dp
    end do
    nodes2(:, 1) = [0.0_dp, 0.0_dp]
    nodes2(:, 2) = [2.0_dp, -1.0_dp]
    nodes2(:, 3) = [3.0_dp, 2.0_dp]
    ! B(s) = [s(4 - s), 2s(2s - 1)]
    expected2(1, :) = s_vals2 * (4.0_dp - s_vals2)
    expected2(2, :) = 2.0_dp * s_vals2 * (2.0_dp * s_vals2 - 1.0_dp)
    call evaluate_multi( &
         3, 2, nodes2, 65, s_vals2, evaluated2)
    case_success = all(evaluated2 == expected2)
    call print_status(name, case_id, case_success, success)

  end subroutine test_evaluate_multi

  subroutine test_specialize_curve(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2)
    real(c_double) :: new_nodes1(2, 2), expected1(2, 2)
    real(c_double) :: nodes2(2, 3), left_nodes(2, 3), right_nodes(2, 3)
    real(c_double) :: new_nodes2(2, 3)
    real(c_double) :: nodes3(2, 4), new_nodes3(2, 4), expected3(2, 4)
    real(c_double) :: nodes4(2, 5), new_nodes4(2, 5), expected4(2, 5)
    integer :: case_id
    character(16) :: name

    case_id = 1
    name = "specialize_curve"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = [0.0_dp, 0.0_dp]
    nodes1(:, 2) = [1.0_dp, 1.0_dp]
    expected1(:, 1) = [0.25_dp, 0.25_dp]
    expected1(:, 2) = [0.75_dp, 0.75_dp]
    call specialize_curve( &
         2, 2, nodes1, 0.25_dp, 0.75_dp, new_nodes1)
    case_success = all(new_nodes1 == expected1)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic curve, after subdivision.
    nodes2(:, 1) = [0.0_dp, 1.0_dp]
    nodes2(:, 2) = [1.0_dp, 6.0_dp]
    nodes2(:, 3) = [3.0_dp, 5.0_dp]
    call subdivide_nodes( &
         3, 2, nodes2, left_nodes, right_nodes)
    call specialize_curve( &
         3, 2, nodes2, 0.0_dp, 0.5_dp, new_nodes2)
    case_success = all(new_nodes2 == left_nodes)
    ! Do a "second" check for the right-hand nodes.
    call specialize_curve( &
         3, 2, nodes2, 0.5_dp, 1.0_dp, new_nodes2)
    case_success = ( &
         case_success .AND. &
         all(new_nodes2 == right_nodes))
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Cubic curve.
    nodes3(:, 1) = [0.0_dp, 0.0_dp]
    nodes3(:, 2) = [1.0_dp, -1.0_dp]
    nodes3(:, 3) = [1.0_dp, -2.0_dp]
    nodes3(:, 4) = [3.0_dp, 2.0_dp]
    expected3(:, 1) = [171.0_dp, -187.0_dp] / 512.0_dp
    expected3(:, 2) = [375.0_dp, -423.0_dp] / 512.0_dp
    expected3(:, 3) = [499.0_dp, -579.0_dp] / 512.0_dp
    expected3(:, 4) = [735.0_dp, -335.0_dp] / 512.0_dp
    call specialize_curve( &
         4, 2, nodes3, 0.125_dp, 0.625_dp, new_nodes3)
    case_success = all(new_nodes3 == expected3)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Quartic curve.
    nodes4(:, 1) = [0.0_dp, 5.0_dp]
    nodes4(:, 2) = [1.0_dp, 6.0_dp]
    nodes4(:, 3) = [1.0_dp, 7.0_dp]
    nodes4(:, 4) = [3.0_dp, 6.0_dp]
    nodes4(:, 5) = [3.0_dp, 7.0_dp]
    expected4(:, 1) = [1.5625_dp, 6.375_dp]
    expected4(:, 2) = [1.78125_dp, 6.4375_dp]
    expected4(:, 3) = [2.015625_dp, 6.46875_dp]
    expected4(:, 4) = [2.2578125_dp, 6.484375_dp]
    expected4(:, 5) = [2.47265625_dp, 6.5234375_dp]
    call specialize_curve( &
         5, 2, nodes4, 0.5_dp, 0.75_dp, new_nodes4)

    case_success = all(new_nodes4 == expected4)
    call print_status(name, case_id, case_success, success)

  end subroutine test_specialize_curve

  subroutine test_evaluate_hodograph(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: s_vals(5), s_val
    real(c_double) :: hodograph(2, 1), expected(2, 1)
    real(c_double) :: nodes1(2, 2), nodes2(2, 3), nodes3(2, 4)
    integer :: i
    integer :: case_id
    character(18) :: name

    case_id = 1
    name = "evaluate_hodograph"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = [0.0_dp, 0.0_dp]
    nodes1(:, 2) = [1.0_dp, 1.0_dp]
    expected(:, 1) = nodes1(:, 2) - nodes1(:, 1)
    s_vals(:2) = [0.25_dp, 0.75_dp]
    case_success = .TRUE.
    do i = 1, 2
       call evaluate_hodograph( &
            s_vals(i), 2, 2, nodes1, hodograph)
       case_success = ( &
            case_success .AND. &
            all(hodograph == expected))
    end do
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic curve.
    nodes2(:, 1) = [0.0_dp, 0.0_dp]
    nodes2(:, 2) = [0.5_dp, 1.0_dp]
    nodes2(:, 3) = [1.25_dp, 0.25_dp]
    !  B(s) = [s(s + 4)/4, s(8 - 7s)/4]
    ! B'(s) = [(2 + s)/2, (4 - 7s)/2]
    s_vals = [0.0_dp, 0.25_dp, 0.5_dp, 0.625_dp, 0.875_dp]
    case_success = .TRUE.
    do i = 1, 5
       s_val = s_vals(i)
       call evaluate_hodograph( &
            s_val, 3, 2, nodes2, hodograph)
       expected(:, 1) = [ &
            (2.0_dp + s_val) / 2.0_dp, &
            (4.0_dp - 7.0_dp * s_val) / 2.0_dp]
       case_success = ( &
            case_success .AND. &
            all(hodograph == expected))
    end do
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Cubic curve.
    nodes3(:, 1) = [0.0_dp, 0.0_dp]
    nodes3(:, 2) = [0.25_dp, 1.0_dp]
    nodes3(:, 3) = [0.75_dp, 0.5_dp]
    nodes3(:, 4) = [1.25_dp, 1.0_dp]
    !  B(s) = [s(3 + 3s - s^2)/4, s(5s^2 - 9s + 6)/2]
    ! B'(s) = [3(1 + 2s - s^2)/4, 3(5s^2 - 6s + 2)/2]
    s_vals = [0.125_dp, 0.5_dp, 0.75_dp, 1.0_dp, 1.125_dp]
    case_success = .TRUE.
    do i = 1, 5
       s_val = s_vals(i)
       call evaluate_hodograph( &
            s_val, 4, 2, nodes3, hodograph)
       expected(1, 1) = ( &
            3.0_dp * (1.0_dp + 2.0_dp * s_val - s_val * s_val) / 4.0_dp)
       expected(2, 1) = ( &
            1.5_dp * (5.0_dp * s_val * s_val - 6.0_dp * s_val + 2.0_dp))
       case_success = ( &
            case_success .AND. &
            all(hodograph == expected))
    end do
    call print_status(name, case_id, case_success, success)

  end subroutine test_evaluate_hodograph

  subroutine test_subdivide_nodes(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2), l_expected1(2, 2), r_expected1(2, 2)
    real(c_double) :: left_nodes1(2, 2), right_nodes1(2, 2)
    real(c_double) :: nodes2(2, 3), l_expected2(2, 3), r_expected2(2, 3)
    real(c_double) :: left_nodes2(2, 3), right_nodes2(2, 3)
    real(c_double) :: nodes3(2, 4), l_expected3(2, 4), r_expected3(2, 4)
    real(c_double) :: left_nodes3(2, 4), right_nodes3(2, 4)
    integer :: case_id
    character(23) :: name

    case_id = 1
    name = "subdivide_nodes (Curve)"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = [0.0_dp, 1.0_dp]
    nodes1(:, 2) = [4.0_dp, 6.0_dp]
    l_expected1(:, 1) = [0.0_dp, 1.0_dp]
    l_expected1(:, 2) = [2.0_dp, 3.5_dp]
    r_expected1(:, 1) = [2.0_dp, 3.5_dp]
    r_expected1(:, 2) = [4.0_dp, 6.0_dp]
    call subdivide_nodes( &
         2, 2, nodes1, left_nodes1, right_nodes1)

    case_success = ( &
         all(left_nodes1 == l_expected1) .AND. &
         all(right_nodes1 == r_expected1))
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Evaluate subdivided parts of a line.
    call subdivide_points_check( &
         2, 2, 909203, 10489, case_success)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Quadratic curve.
    nodes2(:, 1) = [0.0_dp, 1.0_dp]
    nodes2(:, 2) = [4.0_dp, 6.0_dp]
    nodes2(:, 3) = [7.0_dp, 3.0_dp]
    l_expected2(:, 1) = [0.0_dp, 1.0_dp]
    l_expected2(:, 2) = [2.0_dp, 3.5_dp]
    l_expected2(:, 3) = [3.75_dp, 4.0_dp]
    r_expected2(:, 1) = [3.75_dp, 4.0_dp]
    r_expected2(:, 2) = [5.5_dp, 4.5_dp]
    r_expected2(:, 3) = [7.0_dp, 3.0_dp]
    call subdivide_nodes( &
         3, 2, nodes2, left_nodes2, right_nodes2)

    case_success = ( &
         all(left_nodes2 == l_expected2) .AND. &
         all(right_nodes2 == r_expected2))
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Evaluate subdivided parts of a quadratic curve.
    call subdivide_points_check( &
         3, 2, 10920388, 7756, case_success)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Cubic curve.
    nodes3(:, 1) = [0.0_dp, 1.0_dp]
    nodes3(:, 2) = [4.0_dp, 6.0_dp]
    nodes3(:, 3) = [7.0_dp, 3.0_dp]
    nodes3(:, 4) = [6.0_dp, 5.0_dp]
    l_expected3(:, 1) = [0.0_dp, 1.0_dp]
    l_expected3(:, 2) = [2.0_dp, 3.5_dp]
    l_expected3(:, 3) = [3.75_dp, 4.0_dp]
    l_expected3(:, 4) = [4.875_dp, 4.125_dp]
    r_expected3(:, 1) = [4.875_dp, 4.125_dp]
    r_expected3(:, 2) = [6.0_dp, 4.25_dp]
    r_expected3(:, 3) = [6.5_dp, 4.0_dp]
    r_expected3(:, 4) = [6.0_dp, 5.0_dp]
    call subdivide_nodes( &
         4, 2, nodes3, left_nodes3, right_nodes3)

    case_success = ( &
         all(left_nodes3 == l_expected3) .AND. &
         all(right_nodes3 == r_expected3))
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Evaluate subdivided parts of a cubic curve.
    call subdivide_points_check( &
         4, 2, 33981020, 81902302, case_success)
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Quartic curve (not hardcoded).
    call subdivide_points_check( &
         5, 2, 1190, 881, case_success)
    call print_status(name, case_id, case_success, success)

  end subroutine test_subdivide_nodes

  subroutine subdivide_points_check( &
       num_nodes, dimension_, multiplier, modulus, success)
    integer, intent(in) :: num_nodes, dimension_, multiplier, modulus
    logical, intent(out) :: success
    ! Variables outside of signature.
    real(c_double) :: nodes(dimension_, num_nodes)
    real(c_double) :: left(dimension_, num_nodes)
    real(c_double) :: right(dimension_, num_nodes)
    real(c_double) :: unit_interval(33), left_s(33), right_s(33)
    real(c_double) :: evaluated1(dimension_, 33)
    real(c_double) :: evaluated2(dimension_, 33)
    integer :: i

    success = .TRUE.

    call get_random_nodes(nodes, multiplier, modulus, num_bits=8)
    call subdivide_nodes( &
         num_nodes, dimension_, nodes, left, right)

    do i = 1, 33
       unit_interval(i) = (i - 1) / 32.0_dp
       left_s(i) = (i - 1) / 64.0_dp
       right_s(i) = (i + 31) / 64.0_dp
    end do

    call evaluate_multi( &
         num_nodes, dimension_, nodes, 33, left_s, evaluated1)
    call evaluate_multi( &
         num_nodes, dimension_, left, 33, unit_interval, evaluated2)
    success = ( &
         success .AND. &
         all(evaluated1 == evaluated2))

    call evaluate_multi( &
         num_nodes, dimension_, nodes, 33, right_s, evaluated1)
    call evaluate_multi( &
         num_nodes, dimension_, right, 33, unit_interval, evaluated2)
    success = ( &
         success .AND. &
         all(evaluated1 == evaluated2))

  end subroutine subdivide_points_check

  subroutine test_newton_refine(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes(3, 4)
    real(c_double) :: point(3, 1)
    real(c_double) :: new_s
    integer :: case_id
    character(13) :: name

    case_id = 1
    name = "newton_refine"

    ! CASE 1: Cubic in 3D.
    nodes(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
    nodes(:, 2) = [1.0_dp, -1.0_dp, 1.0_dp]
    nodes(:, 3) = [3.0_dp, 2.0_dp, 2.0_dp]
    nodes(:, 4) = [2.0_dp, 2.0_dp, 4.0_dp]
    ! curve(1/2) = p
    point(:, 1) = [1.75_dp, 0.625_dp, 1.625_dp]
    call newton_refine( &
         4, 3, nodes, point, 0.25_dp, new_s)
    case_success = (110.0_dp * new_s == 57.0_dp)
    call print_status(name, case_id, case_success, success)

  end subroutine test_newton_refine

  subroutine test_locate_point(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(3, 3), point1(3, 1)
    real(c_double) :: nodes2(2, 3), point2(2, 1)
    real(c_double) :: nodes3(2, 4), point3(2, 1)
    real(c_double) :: s_val
    integer :: case_id
    character(20) :: name

    case_id = 1
    name = "locate_point (Curve)"

    ! CASE 1: Quadratic in 3D (with match).
    nodes1(:, 1) = 0
    nodes1(:, 2) = [3.0_dp, 0.0_dp, -1.0_dp]
    nodes1(:, 3) = [1.0_dp, 1.0_dp, 3.0_dp]
    point1(:, 1) = [43.0_dp, 1.0_dp, -11.0_dp] / 64.0_dp
    call locate_point( &
         3, 3, nodes1, point1, s_val)

    case_success = (s_val == 0.125_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic in 2D (without match).
    nodes2(:, 1) = 0
    nodes2(:, 2) = [0.5_dp, 1.0_dp]
    nodes2(:, 3) = [1.0_dp, 0.0_dp]
    point2(:, 1) = [0.5_dp, 2.0_dp]
    call locate_point( &
         3, 2, nodes2, point2, s_val)

    case_success = (s_val == LOCATE_MISS)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Self-intersecting cubic in 3D (will cause failure).
    nodes3(:, 1) = [0.0_dp, 2.0_dp]
    nodes3(:, 2) = [-1.0_dp, 0.0_dp]
    nodes3(:, 3) = [1.0_dp, 1.0_dp]
    nodes3(:, 4) = [-0.75_dp, 1.625_dp]
    point3(:, 1) = [-0.25_dp, 1.375_dp]
    call locate_point( &
         4, 2, nodes3, point3, s_val)

    case_success = (s_val == LOCATE_INVALID)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Newton's method pushes the value slightly to the left of ``0.0``.
    nodes2(:, 1) = 0
    nodes2(:, 2) = [1.0_dp, 1.0_dp]
    nodes2(:, 3) = [2.0_dp, 0.0_dp]
    point2(:, 1) = [0.0_dp, 0.0_dp]
    call locate_point( &
         3, 2, nodes2, point2, s_val)

    case_success = (s_val == 0.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Newton's method pushes the value slightly to the right of
    !         ``1.0`` (uses same nodes as CASE 4).
    point2(:, 1) = [2.0_dp, 0.0_dp]
    call locate_point( &
         3, 2, nodes2, point2, s_val)

    case_success = (s_val == 1.0_dp)
    call print_status(name, case_id, case_success, success)

  end subroutine test_locate_point

  subroutine test_elevate_nodes(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2), elevated1(2, 3), expected1(2, 3)
    real(c_double) :: nodes2(3, 3), elevated2(3, 4), expected2(3, 4)
    integer :: case_id
    character(13) :: name

    case_id = 1
    name = "elevate_nodes"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = 0
    nodes1(:, 2) = [2.0_dp, 4.0_dp]
    expected1(:, 1) = [0.0_dp, 0.0_dp]
    expected1(:, 2) = [1.0_dp, 2.0_dp]
    expected1(:, 3) = [2.0_dp, 4.0_dp]
    call elevate_nodes( &
         2, 2, nodes1, elevated1)
    case_success = all(elevated1 == expected1)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic curve.
    nodes2(:, 1) = [0.0_dp, 0.5_dp, 0.75_dp]
    nodes2(:, 2) = [3.0_dp, 0.5_dp, 3.0_dp]
    nodes2(:, 3) = [6.0_dp, 0.5_dp, 2.25_dp]
    expected2(:, 1) = [0.0_dp, 0.5_dp, 0.75_dp]
    expected2(:, 2) = [2.0_dp, 0.5_dp, 2.25_dp]
    expected2(:, 3) = [4.0_dp, 0.5_dp, 2.75_dp]
    expected2(:, 4) = [6.0_dp, 0.5_dp, 2.25_dp]
    call elevate_nodes( &
         3, 3, nodes2, elevated2)
    case_success = all(elevated2 == expected2)
    call print_status(name, case_id, case_success, success)

  end subroutine test_elevate_nodes

  subroutine test_get_curvature(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2), nodes2(2, 3)
    real(c_double) :: s_val, tangent_vec(2, 1), curvature
    integer :: case_id
    character(13) :: name

    case_id = 1
    name = "get_curvature"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = 0
    nodes1(:, 2) = 1
    s_val = 0.5_dp
    call evaluate_hodograph( &
         s_val, 2, 2, nodes1, tangent_vec)
    call get_curvature( &
         2, nodes1, tangent_vec, s_val, curvature)
    case_success = (curvature == 0.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Linear curve, degree-elevated to quadratic.
    nodes2(:, 1) = 0
    nodes2(:, 2) = 0.5_dp
    nodes2(:, 3) = 1
    s_val = 0.25_dp
    call evaluate_hodograph( &
         s_val, 3, 2, nodes2, tangent_vec)
    call get_curvature( &
         3, nodes2, tangent_vec, s_val, curvature)
    case_success = (curvature == 0.0_dp)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Quadratic curve.
    nodes2(:, 1) = 0
    nodes2(:, 2) = [0.5_dp, 1.0_dp]
    nodes2(:, 3) = [1.0_dp, 0.0_dp]
    s_val = 0.5_dp
    call evaluate_hodograph( &
         s_val, 3, 2, nodes2, tangent_vec)
    call get_curvature( &
         3, nodes2, tangent_vec, s_val, curvature)
    case_success = (curvature == -4.0_dp)
    call print_status(name, case_id, case_success, success)

  end subroutine test_get_curvature

  subroutine test_reduce_pseudo_inverse(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    logical(c_bool) :: not_implemented
    real(c_double) :: nodes1(2, 2), reduced1(2, 1), expected1(2, 1)
    real(c_double) :: nodes2(2, 3), reduced2(2, 2), expected2(2, 2)
    real(c_double) :: nodes3(3, 4), reduced3(3, 3), expected3(3, 3)
    real(c_double) :: nodes4(5, 1), reduced4(4, 1), expected4(4, 1)
    real(c_double) :: reduced5(4, 4), expected5(4, 4)
    real(c_double) :: nodes6(2, 6), reduced6(2, 5)
    integer :: case_id
    character(21) :: name

    case_id = 1
    name = "reduce_pseudo_inverse"

    ! CASE 1: Reduce from line to point/constant.
    nodes1(:, 1) = [-2.0_dp, 1.0_dp]
    nodes1(:, 2) = nodes1(:, 1)
    expected1(:, 1) = nodes1(:, 1)
    call reduce_pseudo_inverse( &
         2, 2, nodes1, reduced1, not_implemented)
    case_success = (all(reduced1 == expected1) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Reduce a degree-elevated line (as quadratic).
    nodes2(:, 1) = 0
    nodes2(:, 2) = [1.0_dp, 2.0_dp]
    nodes2(:, 3) = [2.0_dp, 4.0_dp]
    expected2(:, 1) = 0
    expected2(:, 2) = [2.0_dp, 4.0_dp]
    call reduce_pseudo_inverse( &
         3, 2, nodes2, reduced2, not_implemented)
    case_success = (all(reduced2 == expected2) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Test that degree 2->1 reduction is actually a pseudo-inverse.
    call pseudo_inverse_helper(2, reduced2, expected2, not_implemented)
    case_success = (all(reduced2 == expected2) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Reduce a quadratic that is not a deg.-elevated line.
    nodes2(:, 1) = 0
    nodes2(:, 2) = [1.0_dp, 1.5_dp]
    nodes2(:, 3) = [2.0_dp, 0.0_dp]
    expected2(:, 1) = [0.0_dp, 0.5_dp]
    expected2(:, 2) = [2.0_dp, 0.5_dp]
    call reduce_pseudo_inverse( &
         3, 2, nodes2, reduced2, not_implemented)
    case_success = (all(reduced2 == expected2) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Reduce a degree-elevated quadratic (as cubic).
    nodes3(:, 1) = [0.0_dp, 0.5_dp, 0.75_dp]
    nodes3(:, 2) = [2.0_dp, 0.5_dp, 2.25_dp]
    nodes3(:, 3) = [4.0_dp, 0.5_dp, 2.75_dp]
    nodes3(:, 4) = [6.0_dp, 0.5_dp, 2.25_dp]
    expected3(:, 1) = nodes3(:, 1)
    expected3(:, 2) = [3.0_dp, 0.5_dp, 3.0_dp]
    expected3(:, 3) = nodes3(:, 4)
    call reduce_pseudo_inverse( &
         4, 3, nodes3, reduced3, not_implemented)
    case_success = (all(reduced3 == expected3) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Test that degree 3->2 reduction is actually a pseudo-inverse.
    call pseudo_inverse_helper(3, reduced3, expected3, not_implemented)
    case_success = ( &
         maxval(abs(reduced3 - expected3)) < MACHINE_EPS .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Reduce a degree-elevated cubic (as quartic).
    nodes4(:, 1) = [0.0_dp, 0.75_dp, 2.0_dp, 2.75_dp, 2.0_dp]
    expected4(:, 1) = [0.0_dp, 1.0_dp, 3.0_dp, 2.0_dp]
    call reduce_pseudo_inverse( &
         5, 1, nodes4, reduced4, not_implemented)
    case_success = (all(reduced4 == expected4) .AND. .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 8: Test that degree 4->3 reduction is actually a pseudo-inverse.
    call pseudo_inverse_helper(4, reduced5, expected5, not_implemented)
    case_success = ( &
         maxval(abs(reduced3 - expected3)) < MACHINE_EPS .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 9: Unsupported degree
    call get_random_nodes(nodes6, 9177093, 16718, num_bits=8)
    call reduce_pseudo_inverse( &
         6, 2, nodes6, reduced6, not_implemented)
    case_success = (not_implemented)
    call print_status(name, case_id, case_success, success)

  end subroutine test_reduce_pseudo_inverse

  subroutine pseudo_inverse_helper(degree, result_, id_mat, not_implemented)
    integer, intent(in) :: degree
    real(c_double), intent(out) :: result_(degree, degree)
    real(c_double), intent(out) :: id_mat(degree, degree)
    logical(c_bool), intent(out) :: not_implemented
    ! Variables outside of signature.
    real(c_double) :: nodes(degree + 1, degree + 1)
    real(c_double) :: reduction_mat(degree + 1, degree)
    real(c_double) :: elevation_mat(degree, degree + 1)

    nodes = get_id_mat(degree + 1)
    call reduce_pseudo_inverse( &
         degree + 1, degree + 1, nodes, reduction_mat, not_implemented)
    ! NOTE: We assume none of the callers will use inputs that force
    !       ``not_implemeted`` to be ``.FALSE.``.

    id_mat = get_id_mat(degree)
    call elevate_nodes( &
         degree, degree, id_mat, elevation_mat)

    result_ = matmul(elevation_mat, reduction_mat)

  end subroutine pseudo_inverse_helper

  subroutine test_full_reduce(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(1, 2), expected1(1, 2), reduced1(1, 2)
    real(c_double) :: nodes2(2, 4), expected2(2, 4), reduced2(2, 4)
    real(c_double) :: nodes3(2, 6), reduced3(2, 6)
    real(c_double) :: nodes4(2, 1), expected4(2, 1), reduced4(2, 1)
    real(c_double) :: nodes5(2, 5), expected5(2, 5), reduced5(2, 5)
    integer(c_int) :: num_reduced_nodes
    logical(c_bool) :: not_implemented
    integer :: case_id
    character(11) :: name

    case_id = 1
    name = "full_reduce"

    ! CASE 1: Linear curve, one reduction.
    nodes1 = 5.5_dp
    expected1(:, 1) = 5.5_dp
    call full_reduce( &
         2, 1, nodes1, num_reduced_nodes, &
         reduced1, not_implemented)
    case_success = ( &
         num_reduced_nodes == 1 .AND. &
         all(reduced1(:, :1) == expected1(:, :1)) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Cubic curve, one reduction.
    nodes2(:, 1) = 0
    nodes2(:, 2) = [2.0_dp, 4.0_dp]
    nodes2(:, 3) = [4.0_dp, 6.0_dp]
    nodes2(:, 4) = [6.0_dp, 6.0_dp]
    expected2(:, 1) = 0
    expected2(:, 2) = [3.0_dp, 6.0_dp]
    expected2(:, 3) = [6.0_dp, 6.0_dp]
    call full_reduce( &
         4, 2, nodes2, num_reduced_nodes, &
         reduced2, not_implemented)
    case_success = ( &
         num_reduced_nodes == 3 .AND. &
         all(reduced2(:, :3) == expected2(:, :3)) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Cubic curve, two reductions.
    nodes2(:, 1) = [0.0_dp, 4.0_dp]
    nodes2(:, 2) = [1.0_dp, 4.5_dp]
    nodes2(:, 3) = [2.0_dp, 5.0_dp]
    nodes2(:, 4) = [3.0_dp, 5.5_dp]
    expected2(:, 1) = [0.0_dp, 4.0_dp]
    expected2(:, 2) = [3.0_dp, 5.5_dp]
    call full_reduce( &
         4, 2, nodes2, num_reduced_nodes, &
         reduced2, not_implemented)
    case_success = ( &
         num_reduced_nodes == 2 .AND. &
         all(reduced2(:, :2) == expected2(:, :2)) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Cubic curve, no reductions.
    nodes2(:, 1) = [0.0_dp, 2.0_dp]
    nodes2(:, 2) = [-1.0_dp, 0.0_dp]
    nodes2(:, 3) = 1
    nodes2(:, 4) = [-0.75_dp, 1.625_dp]
    expected2 = nodes2
    call full_reduce( &
         4, 2, nodes2, num_reduced_nodes, &
         reduced2, not_implemented)
    case_success = ( &
         num_reduced_nodes == 4 .AND. &
         all(reduced2(:, :4) == expected2(:, :4)) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Unsupported degree (quintic, degree 5).
    nodes3 = 0
    call full_reduce( &
         6, 2, nodes3, num_reduced_nodes, &
         reduced3, not_implemented)
    case_success = (not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Point/constant (can't be any reduction).
    nodes4(:, 1) = [0.0_dp, 1.0_dp]
    expected4 = nodes4
    call full_reduce( &
         1, 2, nodes4, num_reduced_nodes, &
         reduced4, not_implemented)
    case_success = ( &
         num_reduced_nodes == 1 .AND. &
         all(reduced4 == expected4) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Quartic curve, one reduction.
    expected5(:, 1) = [0.0_dp, 0.0_dp]
    expected5(:, 2) = [2.0_dp, 4.0_dp]
    expected5(:, 3) = [4.0_dp, 6.5_dp]
    expected5(:, 4) = [18.0_dp, 1.5_dp]
    nodes5(:, 1) = expected5(:, 1)
    nodes5(:, 2) = (expected5(:, 1) + 3 * expected5(:, 2)) / 4.0_dp
    nodes5(:, 3) = (expected5(:, 2) + expected5(:, 3)) / 2.0_dp
    nodes5(:, 4) = (3 * expected5(:, 3) + expected5(:, 4)) / 4.0_dp
    nodes5(:, 5) = expected5(:, 4)
    call full_reduce( &
         5, 2, nodes5, num_reduced_nodes, &
         reduced5, not_implemented)
    case_success = ( &
         num_reduced_nodes == 4 .AND. &
         all(reduced5(:, :4) == expected5(:, :4)) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

    ! CASE 8: Quartic curve, no reductions.
    nodes5(:, 1) = [0.0_dp, 0.0_dp]
    nodes5(:, 2) = [1.0_dp, 2.0_dp]
    nodes5(:, 3) = [2.0_dp, 3.5_dp]
    nodes5(:, 4) = [5.0_dp, 3.5_dp]
    nodes5(:, 5) = [8.0_dp, 18.0_dp]
    call full_reduce( &
         5, 2, nodes5, num_reduced_nodes, &
         reduced5, not_implemented)
    case_success = ( &
         num_reduced_nodes == 5 .AND. &
         all(reduced5 == nodes5) .AND. &
         .NOT. not_implemented)
    call print_status(name, case_id, case_success, success)

  end subroutine test_full_reduce

  subroutine test_compute_length(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    real(c_double) :: nodes1(2, 2), nodes2(2, 3), nodes3(2, 4)
    real(c_double) :: length, expected
    integer(c_int) :: error_val
    integer :: case_id
    character(14) :: name

    case_id = 1
    name = "compute_length"

    ! CASE 1: Linear curve.
    nodes1(:, 1) = 0
    nodes1(:, 2) = [3.0_dp, 4.0_dp]
    call compute_length( &
         2, 2, nodes1, length, error_val)
    case_success = (length == 5.0_dp .AND. error_val == 0)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Quadratic curve.
    nodes2(:, 1) = 0
    nodes2(:, 2) = [1.0_dp, 2.0_dp]
    nodes2(:, 3) = [2.0_dp, 0.0_dp]
    call compute_length( &
         3, 2, nodes2, length, error_val)
    ! 2 INT_0^1 SQRT(16 s^2  - 16 s + 5) ds = SQRT(5) + sinh^{-1}(2)/2
    expected = sqrt(5.0_dp) + 0.5_dp * asinh(2.0_dp)
    case_success = ( &
         abs(length - expected) <= spacing(expected) .AND. &
         error_val == 0)
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Cubic curve.
    nodes3(:, 1) = 0
    nodes3(:, 2) = [1.0_dp, 2.0_dp]
    nodes3(:, 3) = [2.0_dp, 0.0_dp]
    nodes3(:, 4) = [3.5_dp, 0.0_dp]
    call compute_length( &
         4, 2, nodes3, length, error_val)
    ! x(s) = s (s^2 + 6) / 2
    ! y(s) = 6 s (s - 1)^2
    ! x'(s)^2 + y'(s)^2 = (9/4)(145s^4 - 384s^3 + 356s^2 - 128s + 20)
    ! NOTE: This "literal" has more precision than the type provides.
    expected = 4.0916195514424730685_dp
    case_success = ( &
         abs(length - expected) <= spacing(expected) .AND. &
         error_val == 0)
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Degree-zero "curve".
    call compute_length( &
         1, 2, nodes1(:, :1), length, error_val)
    case_success = (length == 0.0_dp .AND. error_val == 0)
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Empty curve.
    call compute_length( &
         0, 2, nodes1(:, :0), length, error_val)
    case_success = error_val == -1
    call print_status(name, case_id, case_success, success)

  end subroutine test_compute_length

  subroutine test_curves_equal(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer :: case_id
    type(CurveData) :: &
         curve1, curve2, curve3, curve4, curve5, curve6, curve7
    character(12) :: name

    case_id = 1
    name = "curves_equal"

    ! curve1 is the basis for all comparisons.
    curve1%start = 0.25_dp
    curve1%end_ = 0.75_dp
    allocate(curve1%nodes(2, 2))
    curve1%nodes(:, 1) = [0.5_dp, 2.5_dp]
    curve1%nodes(:, 2) = [7.0_dp, -2.0_dp]

    ! CASE 1: Compare curve1 to itself.
    case_success = curves_equal(curve1, curve1)
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Compare curve1 to uninitialized nodes.
    curve2%start = curve1%start
    curve2%end_ = curve1%end_
    case_success = ( &
         .NOT. allocated(curve2%nodes) .AND. &
         .NOT. curves_equal(curve1, curve2) .AND. &
         .NOT. curves_equal(curve2, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 3: Compare curve1 to an identical, but different curve.
    curve3%start = curve1%start
    curve3%end_ = curve1%end_
    curve3%nodes = curve1%nodes
    case_success = ( &
         curves_equal(curve1, curve3) .AND. curves_equal(curve3, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 4: Compare curve1 to a curve that differs only in ``start``.
    curve4%start = curve1%start - 0.25_dp
    curve4%end_ = curve1%end_
    curve4%nodes = curve1%nodes
    case_success = ( &
         .NOT. curves_equal(curve1, curve4) .AND. &
         .NOT. curves_equal(curve4, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 5: Compare curve1 to a curve that differs only in ``end_``.
    curve5%start = curve1%start
    curve5%end_ = curve1%end_ + 0.25_dp
    curve5%nodes = curve1%nodes
    case_success = ( &
         .NOT. curves_equal(curve1, curve5) .AND. &
         .NOT. curves_equal(curve5, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 6: Compare curve1 to a curve that differs in node shape
    curve6%start = curve1%start
    curve6%end_ = curve1%end_
    allocate(curve6%nodes(2, 1))
    curve6%nodes(:, 1) = curve1%nodes(:, 1)
    case_success = ( &
         .NOT. curves_equal(curve1, curve6) .AND. &
         .NOT. curves_equal(curve6, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 7: Compare curve1 to a curve that differs in
    !         ``nodes`` value but not shape.
    curve7%start = curve1%start
    curve7%end_ = curve1%end_
    curve7%nodes = curve1%nodes + 2.5_dp
    case_success = ( &
         .NOT. curves_equal(curve1, curve7) .AND. &
         .NOT. curves_equal(curve7, curve1))
    call print_status(name, case_id, case_success, success)

    ! CASE 8: Compare curve with uninitialized nodes to **itself**.
    case_success = ( &
         .NOT. allocated(curve2%nodes) .AND. &
         curves_equal(curve2, curve2))
    call print_status(name, case_id, case_success, success)

  end subroutine test_curves_equal

  subroutine test_subdivide_curve(success)
    logical(c_bool), intent(inout) :: success
    ! Variables outside of signature.
    logical :: case_success
    integer :: case_id
    type(CurveData) :: curve_data, left1, left2, right
    real(c_double) :: left_nodes(2, 4), right_nodes(2, 4)
    character(15) :: name

    case_id = 1
    name = "subdivide_curve"

    ! CASE 1: Curve has not yet been subdivided.
    allocate(curve_data%nodes(2, 4))
    curve_data%nodes(:, 1) = [0.0_dp, 1.0_dp]
    curve_data%nodes(:, 2) = [1.0_dp, 2.0_dp]
    curve_data%nodes(:, 3) = [1.0_dp, 2.0_dp]
    curve_data%nodes(:, 4) = [3.0_dp, 0.0_dp]

    left_nodes(:, 1) = [0.0_dp, 1.0_dp]
    left_nodes(:, 2) = [0.5_dp, 1.5_dp]
    left_nodes(:, 3) = [0.75_dp, 1.75_dp]
    left_nodes(:, 4) = [1.125_dp, 1.625_dp]
    right_nodes(:, 1) = [1.125_dp, 1.625_dp]
    right_nodes(:, 2) = [1.5_dp, 1.5_dp]
    right_nodes(:, 3) = [2.0_dp, 1.0_dp]
    right_nodes(:, 4) = [3.0_dp, 0.0_dp]

    call subdivide_curve(curve_data, left1, right)
    case_success = ( &
         left1%start == 0.0_dp .AND. &
         left1%end_ == 0.5_dp .AND. &
         all(left1%nodes == left_nodes) .AND. &
         right%start == 0.5_dp .AND. &
         right%end_ == 1.0_dp .AND. &
         all(right%nodes == right_nodes))
    call print_status(name, case_id, case_success, success)

    ! CASE 2: Curve has **already** been subdivided.
    left_nodes(:, 1) = [0.0_dp, 1.0_dp]
    left_nodes(:, 2) = [0.25_dp, 1.25_dp]
    left_nodes(:, 3) = [0.4375_dp, 1.4375_dp]
    left_nodes(:, 4) = [0.609375_dp, 1.546875_dp]
    right_nodes(:, 1) = [0.609375_dp, 1.546875_dp]
    right_nodes(:, 2) = [0.78125_dp, 1.65625_dp]
    right_nodes(:, 3) = [0.9375_dp, 1.6875_dp]
    right_nodes(:, 4) = [1.125_dp, 1.625_dp]

    call subdivide_curve(left1, left2, right)
    case_success = ( &
         left2%start == 0.0_dp .AND. &
         left2%end_ == 0.25_dp .AND. &
         all(left2%nodes == left_nodes) .AND. &
         right%start == 0.25_dp .AND. &
         right%end_ == 0.5_dp .AND. &
         all(right%nodes == right_nodes))
    call print_status(name, case_id, case_success, success)

  end subroutine test_subdivide_curve

end module test_curve

max_stars_repo_name: dibir-magomedsaygitov/bezier

max_issues_repo_name: dibir-magomedsaygitov/bezier

max_forks_repo_name: dibir-magomedsaygitov/bezier


text:
Koko has nothing to say at the moment. Sorry.

20090915 00:10:14 nbsp When you took the link to Tasha off you almost Orphaned Pages orphaned the page for her. Users/JasonAller

20091226 14:37:09 nbsp Clue #5:

I go up and down,
yet never seem to move.
To find the next clue,
look near the hand groove. Users/MikeMcClain

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
      PROGRAM SOMP
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION AA(5, 5)
      DIMENSION BB(5,5)
C

      CALL SOMTAB(AA,BB)
C
      END
C
C     On doit obtenir une region MAY pour A
C     et une region MUST pour B, et ces tags
C     doivent etre conserve's pour les summary regions
C
      SUBROUTINE SOMTAB(A,B)
C
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION A(5,5)
      DIMENSION B(5,5)
C
      B(1,1) = 3
      M = A(1,1) + A(3,3)
C
      END


text:
program pi_main

  use timer
  use omp_lib

  ! Local variables
  integer, parameter :: num_steps = 100000000 ! number of steps over which to estimate pi
  real(kind=8) :: step                        ! the step size
  integer :: ii                               ! genereric counter
  real(kind=8) :: x                           ! intermediate value
  real(kind=8) :: pi = 0.0_8                  ! overall estimate
  real(kind=8), allocatable :: sum(:)         ! variable to store partial sum
  real(kind=8) :: start, end                  ! timers
  integer :: nthreads                         ! number of OpenMP threads
  integer :: tid                              ! thread id

  real(kind=8), parameter :: PI_8 = 4.0_8 * atan(1.0_8)

  ! Get number of OpenMP threads
  !$omp parallel
  nthreads = omp_get_num_threads()
  !$omp end parallel

  allocate(sum(nthreads))

  ! step size is dependent upon the number of steps
  step = 1.0_8/num_steps

  ! Start timer
  call wtime(start)

  ! main loop
  !$omp parallel private(x,tid)
  tid = omp_get_thread_num()
  sum(tid+1) = 0.0_8
  !$omp do
  do ii = 1, num_steps
    x = (ii-0.5_8)*step
    sum(tid+1) = sum(tid+1) + (4.0_8/(1.0_8+x*x))
    !$omp flush(sum)
  end do
  !$omp end do
  !$omp end parallel

  ! Total partial sums serially
  do ii = 1, nthreads
    pi = pi + sum(ii)
  end do
  pi = pi * step

  ! Stop timer
  call wtime(end)

  ! Print result
  write(*,"(A)")         "------------------------------------"
  write(*,"(A,F19.16)") "pi is:    ", pi
  write(*,"(A,F19.16)") "error is: ", abs(pi - PI_8)
  write(*,"(A,F10.3)")   "runtime:  ", end-start
  write(*,"(A)")         "------------------------------------"

  deallocate(sum)

end program pi_main

max_stars_repo_name: PourroyJean/openmp-for-cs

max_issues_repo_name: PourroyJean/openmp-for-cs

max_forks_repo_name: PourroyJean/openmp-for-cs


text:
module trxnvars !  Basic variables shared among trxn's subroutines
!***********************************************************************
! Copyright 2012. Los Alamos National Security, LLC.  This material was
! produced under U.S. Government contract DE-AC52-06NA25396 for Los
! Alamos National Laboratory (LANL), which is operated by Los Alamos
! National Security, LLC for the U.S. Department of Energy. The U.S.
! Government has rights to use, reproduce, and distribute this software.
! Neither the U.S. Government nor Los Alamos National Security, LLC or
! persons acting on their behalf, make any warranty, express or implied,
! or assumes any liability for the accuracy, completeness, or usefulness
! of the software, any information pertaining to the software, or
! represents that its use would not infringe privately owned rights.
!
! The software being licensed may be Export Controlled.   It may not be
! distributed or used by individuals or entities prohibited from having
! access to the software package, pursuant to United States export
! control laws and regulations.  An export control review and
! determination must be completed before LANS will provide access to the
! identified Software.
!***********************************************************************

        use comchem, only : NAME_MAX
	integer SPEC_MAX, LINE_MAX		! Parameters for the maximum space to allocate for different fields.  Probably should get rid of them
	character*4 aSPEC_MAX, aNAME_MAX
	parameter ( SPEC_MAX = 100 )
	parameter ( LINE_MAX = 200 )
	parameter ( aSPEC_MAX = '100' )
	parameter ( aNAME_MAX = '40' )
	integer inpttmp					! Holder for input file handle
	integer linen					! Number of the current line
	character(LINE_MAX) line			! The current line being read (scrapped)
	logical debug             			! Whether we are in debug mode
	logical debug_stop				! Whether to halt execution after reading the macro
	integer trxn_flag         			! Whether we are calling tracrxn or trac
	integer numzones, zonemax			! Number of zones and maximum zone number (persistent)
	integer, allocatable :: zonenums(:)		! Numbers corresponding to zones
	character(20), allocatable :: zonenames(:)	! Zone names
	logical molein					! Whether we are using the total mole initial water option
	logical interactive				! Whether we are reading interactively

! lookupvars:  Variables for database lookup subroutine in trxn
        logical flookup                                                 ! Whether lookup has been called
        character(LINE_MAX), dimension(SPEC_MAX) :: options             ! Options
        character(LINE_MAX) database                                    ! Database file name
        integer dat                                                     ! Filehandle for database
        integer nlocomp, nlocompaq                                      ! Number of components lookup needs to add
        character(NAME_MAX), dimension(SPEC_MAX, 3) :: locomp           ! Components to be added by lookup (1 = state, 2 = component, 3 = master)
        integer ndatcomp                                                ! Number of complexes in the database
        character(NAME_MAX), allocatable :: datcomp(:), datmaster(:)    ! Complexes in the database and their master species
        integer ndatcplx, ndatmin, maxdatcplx, maxdatmin                ! Number of complexes and minerals in the database, and the maximum number of reactants and products for each
        real*8, allocatable :: datcden(:)                               ! Molecular weights in the database
        character(NAME_MAX), allocatable :: datcplxmain(:)              ! Complex formulae
        character(NAME_MAX), allocatable :: datcplx(:, :)               ! Species for complexes in the database.  datcplx(i, j) is reaction i, species j
        real*8, allocatable :: datcplxstoic(:, :)                       ! Stoichiometry for complexes in the database.  Same indexing as datcplx
        real*8, allocatable :: datcplxlkeq(:)                           ! log_k for complexes in the database
        real*8, allocatable :: datcplxheq(:)                            ! delta_h for complexes in the database
        real*8, allocatable :: datcplxtemp(:, :)                        ! Temperature dependency coefficients for complexes in the database
        character(NAME_MAX), allocatable :: datminnam(:)                ! Names of minerals in the database
        character(NAME_MAX), allocatable :: datminmain(:)               ! Mineral formulae
        character(NAME_MAX), allocatable :: datmin(:, :)                ! Species of minerals in the database.  Same indexing as datcplx
        real*8, allocatable :: datminstoic(:, :)                        ! Stoichiometry for minerals in the database. Same indexing as datcplx
        real*8, allocatable :: datminlkeq(:)                            ! log_k for minerals in the database
        real*8, allocatable :: datminheq(:)                             ! delta_h for minerals in the database
        real*8, allocatable :: datmintemp(:, :)                         ! Temperature dependency coefficients for minerals in the database

end module trxnvars

max_stars_repo_name: satkarra/FEHM

max_issues_repo_name: satkarra/FEHM

max_forks_repo_name: satkarra/FEHM


text:
Describe Users/KimberlyMorris here.

20110209 15:01:00 nbsp Welcome to the Wiki/Business Owner Welcome to the Wiki!, Kimberly! Sorry about editing over youwas just cleaning up a bit. Users/PeterBoulay

20110209 15:09:34 nbsp KimI moved the img up and to the right to take advantage of the blank space to the right of the headers Users/PeterBoulay

20110209 15:22:13 nbsp Picture looks great from here. :) Users/PeterBoulay


text:
!==============================================================================!
  subroutine Compute_Gradient_Matrix(grid, Boundary)
!------------------------------------------------------------------------------!
!   Calculates gradient matrix.                                                !
!------------------------------------------------------------------------------!
!----------------------------------[Modules]-----------------------------------!
  use all_mod
  use pro_mod
  use Grid_Mod
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  type(Grid_Type) :: grid
  logical         :: Boundary
!-----------------------------------[Locals]-----------------------------------!
  integer :: c, c1, c2, s
  real    :: Dxc1, Dyc1, Dzc1, Dxc2, Dyc2, Dzc2
  real    :: Jac, Ginv(6)
!==============================================================================!

  do c = 1, grid % n_cells
    G(1,c) = 0.0
    G(2,c) = 0.0
    G(3,c) = 0.0
    G(4,c) = 0.0
    G(5,c) = 0.0
    G(6,c) = 0.0
  end do

  do s = 1, grid % n_faces
    c1 = grid % faces_c(1,s)
    c2 = grid % faces_c(2,s)

    Dxc1 = grid % dx(s)
    Dyc1 = grid % dy(s)
    Dzc1 = grid % dz(s)
    Dxc2 = grid % dx(s)
    Dyc2 = grid % dy(s)
    Dzc2 = grid % dz(s)

    ! Take care of material interfaces         ! 2mat
    if( StateMat(material(c1))==FLUID .and. &  ! 2mat
        StateMat(material(c2))==SOLID       &  ! 2mat
        .or.                                &  ! 2mat
        StateMat(material(c1))==SOLID .and. &  ! 2mat
        StateMat(material(c2))==FLUID ) then   ! 2mat
      Dxc1 = grid % xf(s) - grid % xc(c1)      ! 2mat
      Dyc1 = grid % yf(s) - grid % yc(c1)      ! 2mat
      Dzc1 = grid % zf(s) - grid % zc(c1)      ! 2mat
      Dxc2 = grid % xf(s) - grid % xc(c2)      ! 2mat
      Dyc2 = grid % yf(s) - grid % yc(c2)      ! 2mat
      Dzc2 = grid % zf(s) - grid % zc(c2)      ! 2mat
    end if                                     ! 2mat

    ! With boundary cells, velocities, temperatures
    if(Boundary) then
      G(1,c1)=G(1,c1) + Dxc1*Dxc1  ! 1,1
      G(2,c1)=G(2,c1) + Dyc1*Dyc1  ! 2,2
      G(3,c1)=G(3,c1) + Dzc1*Dzc1  ! 3,3
      G(4,c1)=G(4,c1) + Dxc1*Dyc1  ! 1,2  &  2,1
      G(5,c1)=G(5,c1) + Dxc1*Dzc1  ! 1,3  &  3,1
      G(6,c1)=G(6,c1) + Dyc1*Dzc1  ! 2,3  &  3,2
      if(c2  > 0) then  ! this is enough even for parallel
        G(1,c2)=G(1,c2) + Dxc2*Dxc2  ! 1,1
        G(2,c2)=G(2,c2) + Dyc2*Dyc2  ! 2,2
        G(3,c2)=G(3,c2) + Dzc2*Dzc2  ! 3,3
        G(4,c2)=G(4,c2) + Dxc2*Dyc2  ! 1,2  &  2,1
        G(5,c2)=G(5,c2) + Dxc2*Dzc2  ! 1,3  &  3,1
        G(6,c2)=G(6,c2) + Dyc2*Dzc2  ! 2,3  &  3,2
      end if

    ! Without boundary cells => pressure
    else ! Don't use Boundary
      if(c2 > 0 .or. c2 < 0 .and. TypeBC(c2) == BUFFER) then
        G(1,c1)=G(1,c1) + Dxc1*Dxc1  ! 1,1
        G(2,c1)=G(2,c1) + Dyc1*Dyc1  ! 2,2
        G(3,c1)=G(3,c1) + Dzc1*Dzc1  ! 3,3
        G(4,c1)=G(4,c1) + Dxc1*Dyc1  ! 1,2  &  2,1
        G(5,c1)=G(5,c1) + Dxc1*Dzc1  ! 1,3  &  3,1
        G(6,c1)=G(6,c1) + Dyc1*Dzc1  ! 2,3  &  3,2
      end if
      if(c2 > 0) then
        G(1,c2)=G(1,c2) + Dxc2*Dxc2  ! 1,1
        G(2,c2)=G(2,c2) + Dyc2*Dyc2  ! 2,2
        G(3,c2)=G(3,c2) + Dzc2*Dzc2  ! 3,3
        G(4,c2)=G(4,c2) + Dxc2*Dyc2  ! 1,2  &  2,1
        G(5,c2)=G(5,c2) + Dxc2*Dzc2  ! 1,3  &  3,1
        G(6,c2)=G(6,c2) + Dyc2*Dzc2  ! 2,3  &  3,2
      end if
    end if ! Boundary
  end do

  !----------------------------------!
  !   Find the inverse of matrix G   !
  !----------------------------------!
  do c = 1, grid % n_cells
    Jac  =         G(1,c) * G(2,c) * G(3,c)  &
           -       G(1,c) * G(6,c) * G(6,c)  &
           -       G(4,c) * G(4,c) * G(3,c)  &
           + 2.0 * G(4,c) * G(5,c) * G(6,c)  &
           -       G(5,c) * G(5,c) * G(2,c)

    Ginv(1) = +( G(2,c)*G(3,c) - G(6,c)*G(6,c) ) / (Jac+TINY)
    Ginv(2) = +( G(1,c)*G(3,c) - G(5,c)*G(5,c) ) / (Jac+TINY)
    Ginv(3) = +( G(1,c)*G(2,c) - G(4,c)*G(4,c) ) / (Jac+TINY)
    Ginv(4) = -( G(4,c)*G(3,c) - G(5,c)*G(6,c) ) / (Jac+TINY)
    Ginv(5) = +( G(4,c)*G(6,c) - G(5,c)*G(2,c) ) / (Jac+TINY)
    Ginv(6) = -( G(1,c)*G(6,c) - G(4,c)*G(5,c) ) / (Jac+TINY)

    G(1,c) = Ginv(1)
    G(2,c) = Ginv(2)
    G(3,c) = Ginv(3)
    G(4,c) = Ginv(4)
    G(5,c) = Ginv(5)
    G(6,c) = Ginv(6)
  end do

  end subroutine

max_stars_repo_name: mhadziabdic/T-FlowS

max_issues_repo_name: mhadziabdic/T-FlowS

max_forks_repo_name: mhadziabdic/T-FlowS