Добавили папку zve

.gitignore

@@ -12,3 +12,4 @@
 !memcached/
 !gravity/
 !fte/
+!zve/

zve/.gitignore

+*
+write_buildinfo.f95
+!description.txt
+!buildinfo
+!eph/JPLEPH
+!.gitignore
+!*.f
+!*.f95
+!*.for
+!*.cmn
+!*.inc
+!*.sh
+!makefile
+!Makefile
+!eph/
+!mod/
+!math/
+!modules/
+!gencatalog/

zve/Makefile

+PLATFORM:=$(shell uname -s |awk -F- '{print $$1}')
+
+FC=gfortran
+RM=rm
+AR = ar
+ARFLAGS= rc
+
+BUILDDIR = ../build/$(PLATFORM)
+MODDIR = $(BUILDDIR)/mod
+LIBDIR = $(BUILDDIR)/lib
+SHAREDIR = $(BUILDDIR)/share
+OBJDIR = $(BUILDDIR)/obj/solvers/solv
+TSTDIR = $(BUILDDIR)/test
+
+JPL=eph/JPLEPH
+JPLTARG = $(subst eph/,$(SHAREDIR)/,$(JPL))
+
+LDFLAGS  = -lrdjson -lwrjson -lfson -lsolv -lsofa -lmemc_gate -lgravity -lmemcached -lsofa_c -lstdc++ -L$(LIBDIR)
+FCFLAGS  = -MD -Wall -fmax-stack-var-size=3932160 -ffree-line-length-none -O -J$(MODDIR) -cpp
+
+LIBSRC = \
+    $(shell ls eph/*.f*) \
+    $(shell ls math/*.f*) \
+    $(shell ls library/*.f*) \
+    $(shell ls mod/*.f*)
+
+OBJS = \
+	$(subst eph/,$(OBJDIR)/, \
+	$(subst math/,$(OBJDIR)/, \
+	$(subst library/,$(OBJDIR)/, \
+	$(subst mod/,$(OBJDIR)/, \
+		$(patsubst %.f,%.o,$(patsubst %.f95,%.o,$(LIBSRC)))))))
+
+DEPS = \
+	$(subst eph/,$(OBJDIR)/, \
+	$(subst math/,$(OBJDIR)/, \
+	$(subst library/,$(OBJDIR)/, \
+	$(subst mod/,$(OBJDIR)/, \
+	    $(patsubst %.f,%.d,$(patsubst %.f95,%.d,$(LIBSRC)))))))
+
+
+MNMODS = $(wildcard ../mod/*.mod)
+MODSRC = $(wildcard modules/*.f*)
+MODOBJ = $(subst modules/,$(OBJDIR)/,$(patsubst %.f95,%.o,$(MODSRC)))
+MODDEP = $(subst modules/,$(OBJDIR)/,$(patsubst %.f95,%.d,$(MODSRC)))
+MODMOD = $(subst modules/,$(MODDIR)/,$(patsubst %.f95,%.mod,$(MODSRC)))
+
+
+LIBTARGET=$(LIBDIR)/libsolv.a
+
+ifndef DEBUG
+	RELCFLAGS = -O2
+else
+	RELCFLAGS = -g -O0 -DDEBUG=$(DEBUG)
+endif
+
+SOLVERS=gencatalog 
+
+all: dirs lib $(SOLVERS)
+	cp $(JPL) $(JPLTARG)
+
+.PHONY: subdirs $(SOLVERS)
+
+subdirs: $(SOLVERS)
+
+$(SOLVERS):
+	$(MAKE) -C $@
+
+lib: dirs $(LIBTARGET)
+
+dirs:
+	@mkdir -p $(MODDIR) $(LIBDIR) $(TSTDIR) $(SHAREDIR) $(OBJDIR)
+
+$(LIBTARGET) : $(MODOBJ) $(OBJS) $(DEPS) $(LIBSRC) $(MODSRC) $(MNMODS)
+	$(AR) $(ARFLAGS) $(LIBTARGET) $(MODOBJ) $(OBJS)
+
+$(OBJDIR)/%.o: modules/%.f*
+	$(FC)  $(FCFLAGS) $(RELCFLAGS) -c $< -o $@
+
+$(OBJDIR)/%.o: math/%.f*
+	$(FC)  $(FCFLAGS) $(RELCFLAGS) -c $< -o $@
+
+$(OBJDIR)/%.o: mod/%.f*
+	$(FC)  $(FCFLAGS) $(RELCFLAGS) -c $< -o $@
+
+$(OBJDIR)/%.o: library/%.f* $(MNMODS)
+	$(FC)  $(FCFLAGS) $(RELCFLAGS) -c $< -o $@
+
+$(OBJDIR)/%.o: eph/%.f*
+	$(FC)  $(FCFLAGS) $(RELCFLAGS) -c $< -o $@
+
+include $(wildcard $(OBJDIR)/*.d)
+
+
+clean:
+	@rm -rf $(OBJS) $(DEPS) $(LIBTARGET) $(MODMOD) $(MODOBJ) $(MODDEP)
+	-for d in $(SOLVERS); do (cd $$d; $(MAKE) clean ); done
+

zve/buildinfo

+#!/bin/sh
+
+
+if [ ! -f __build__ ]; then
+	echo "0" >__build__
+fi
+cr=`cat __build__` 
+num=`expr $cr + 1`
+echo $num > __build__
+
+echo '!> @brief Вывод информации об исполняемомом файле.
+!> @details Подпрограмма выводит в стандартный поток stderr информацию о:
+!>    - номере версии сборки
+!>    - дате сборки исполяемого файла
+!>    - версии компилятора
+!>    - номере "снимка" git
+!>
+!> @note код обновляется автоматически при компиляции.
+subroutine write_buildinfo()
+    use, intrinsic :: iso_fortran_env
+    write(ERROR_UNIT,'\''(A,I0.0)'\'') " LOG: Version: 3.8.",'$num'
+    write(ERROR_UNIT,*) "LOG: Build: ",__DATE__," ",__TIME__
+    write(ERROR_UNIT,*) "LOG: Compiler: ",__VERSION__
+    write(ERROR_UNIT,*) "LOG: Commit: ","'`git rev-parse HEAD`'"
+end subroutine
+' > write_buildinfo.f95
+

zve/eph/fsizer.f

+      SUBROUTINE FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
+C
+C++++++++++++++++++++++++
+C
+C   Version 1.0 uses the INQUIRE statement to find out the the record length 
+C   of the direct access file before opening it.  This procedure is non-standard, 
+C   but seems to work for VAX machines. 
+C
+C  THE SUBROUTINE ALSO SETS THE VALUES OF  NRECL, NRFILE, AND NAMFIL.
+
+C  *****************************************************************
+C  *****************************************************************
+C
+C  THE PARAMETERS NAMFIL, NRECL, AND NRFILE ARE TO BE SET BY THE USER
+C
+C  *****************************************************************
+
+C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
+
+      CHARACTER*255  NAMFIL
+
+C  *****************************************************************
+
+C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
+C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
+C  (for a VAX, it is probably 1)
+C
+      NRECL=4 
+
+C  *****************************************************************
+
+C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
+
+      NRFILE=12
+
+C  *****************************************************************
+
+C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
+
+      CALL getenv('JPLEPH',NAMFIL)
+
+C  *****************************************************************
+
+C   FIND THE RECORD SIZE USING THE INQUIRE STATEMENT  
+
+
+      IRECSZ=0
+
+      INQUIRE(FILE=NAMFIL,RECL=IRECSZ)
+
+C IF 'INQUIRE' DOES NOT WORK, USUALLY IRECSZ WILL BE LEFT AT 0 
+
+      IF(IRECSZ .LE. 0) write(*,*)
+     . ' INQUIRE STATEMENT PROBABLY DID NOT WORK'
+
+      KSIZE=IRECSZ/NRECL
+
+      RETURN
+
+      END
+C++++++++++++++++++++++++
+C
+      SUBROUTINE FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
+C
+C++++++++++++++++++++++++
+C  THIS SUBROUTINE OPENS THE FILE, 'NAMFIL', WITH A PHONY RECORD LENGTH, READS 
+C  THE FIRST RECORD, AND USES THE INFO TO COMPUTE KSIZE, THE NUMBER OF SINGLE 
+C  PRECISION WORDS IN A RECORD.  
+C
+C  THE SUBROUTINE ALSO SETS THE VALUES OF  NRECL, NRFILE, AND NAMFIL.
+
+      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
+
+      SAVE
+
+      CHARACTER*6 TTL(14,3),CNAM(400)
+      CHARACTER*255 NAMFIL
+
+      DIMENSION SS(3)
+
+      INTEGER IPT(3,13)
+
+C  *****************************************************************
+C  *****************************************************************
+C
+C  THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
+C
+C  *****************************************************************
+
+C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
+C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
+C  (for UNIX, it is probably 4)
+C
+      NRECL= 4
+
+C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE
+
+      NRFILE=12
+
+C  *****************************************************************
+
+C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
+
+      CALL getenv('JPLEPH',NAMFIL)
+
+
+C  *****************************************************************
+C  *****************************************************************
+
+C  **  OPEN THE DIRECT-ACCESS FILE AND GET THE POINTERS IN ORDER TO 
+C  **  DETERMINE THE SIZE OF THE EPHEMERIS RECORD
+
+      MRECL=NRECL*1000
+
+        OPEN(NRFILE,
+     *       FILE=NAMFIL,
+     *       ACCESS='DIRECT',
+     *       FORM='UNFORMATTED',
+     *       RECL=MRECL,
+     *       STATUS='OLD')
+
+      READ(NRFILE,REC=1)TTL,CNAM,SS,NCON,AU,EMRAT,
+     * ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
+
+      CLOSE(NRFILE)
+
+
+
+C  FIND THE NUMBER OF EPHEMERIS COEFFICIENTS FROM THE POINTERS
+
+      KMX = 0
+      KHI = 0
+
+      DO I = 1,13
+         IF (IPT(1,I) .GT. KMX) THEN
+            KMX = IPT(1,I)
+            KHI = I
+         ENDIF
+      ENDDO
+
+      ND = 3
+      IF (KHI .EQ. 12) ND=2
+
+      KSIZE = 2*(IPT(1,KHI)+ND*IPT(2,KHI)*IPT(3,KHI)-1)
+
+      RETURN
+
+      END
+C++++++++++++++++++++++++
+C
+      SUBROUTINE FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
+C
+C++++++++++++++++++++++++
+C
+C  THE SUBROUTINE SETS THE VALUES OF  NRECL, KSIZE, NRFILE, AND NAMFIL.
+
+      SAVE
+
+      CHARACTER*80 NAMFIL
+
+C  *****************************************************************
+C  *****************************************************************
+C
+C  THE PARAMETERS NRECL, NRFILE, AND NAMFIL ARE TO BE SET BY THE USER
+
+C  *****************************************************************
+
+C  NRECL=1 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN S.P. WORDS
+C  NRECL=4 IF "RECL" IN THE OPEN STATEMENT IS THE RECORD LENGTH IN BYTES
+
+      NRECL= 4
+
+C  *****************************************************************
+
+C  NRFILE IS THE INTERNAL UNIT NUMBER USED FOR THE EPHEMERIS FILE (DEFAULT: 12)
+
+      NRFILE=12
+
+C  *****************************************************************
+
+C  NAMFIL IS THE EXTERNAL NAME OF THE BINARY EPHEMERIS FILE
+
+      CALL getenv('JPLEPH',NAMFIL)
+C      NAMFIL='JPLEPH'
+
+C  *****************************************************************
+
+C  KSIZE must be set by the user according to the ephemeris to be read
+
+C  For  de200, set KSIZE to 1652
+C  For  de405, set KSIZE to 2036
+C  For  de406, set KSIZE to 1456
+C  For  de414, set KSIZE to 2036
+
+      KSIZE =  2036
+
+C  *******************************************************************
+
+      RETURN
+
+      END

zve/eph/interp.f

+      SUBROUTINE INTERP(BUF,T,NCF,NCM,NA,IFL,PV)
+C
+C+++++++++++++++++++++++++++++++++
+C
+C     THIS SUBROUTINE DIFFERENTIATES AND INTERPOLATES A
+C     SET OF CHEBYSHEV COEFFICIENTS TO GIVE POSITION AND VELOCITY
+C
+C     CALLING SEQUENCE PARAMETERS:
+C
+C       INPUT:
+C
+C         BUF   1ST LOCATION OF ARRAY OF D.P. CHEBYSHEV COEFFICIENTS OF POSITION
+C
+C           T   T(1) IS DP FRACTIONAL TIME IN INTERVAL COVERED BY
+C               COEFFICIENTS AT WHICH INTERPOLATION IS WANTED
+C               (0 .LE. T(1) .LE. 1).  T(2) IS DP LENGTH OF WHOLE
+C               INTERVAL IN INPUT TIME UNITS.
+C
+C         NCF   # OF COEFFICIENTS PER COMPONENT
+C
+C         NCM   # OF COMPONENTS PER SET OF COEFFICIENTS
+C
+C          NA   # OF SETS OF COEFFICIENTS IN FULL ARRAY
+C               (I.E., # OF SUB-INTERVALS IN FULL INTERVAL)
+C
+C          IFL  INTEGER FLAG: =1 FOR POSITIONS ONLY
+C                             =2 FOR POS AND VEL
+C
+C
+C       OUTPUT:
+C
+C         PV   INTERPOLATED QUANTITIES REQUESTED.  DIMENSION
+C               EXPECTED IS PV(NCM,IFL), DP.
+C
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+C
+      SAVE
+C
+      REAL*8 BUF(NCF,NCM,*),T(2),PV(NCM,*),PC(18),VC(18)
+
+C
+      DATA NP/2/
+      DATA NV/3/
+      DATA TWOT/0.D0/
+      DATA PC(1),PC(2)/1.D0,0.D0/
+      DATA VC(2)/1.D0/
+C
+C       ENTRY POINT. GET CORRECT SUB-INTERVAL NUMBER FOR THIS SET
+C       OF COEFFICIENTS AND THEN GET NORMALIZED CHEBYSHEV TIME
+C       WITHIN THAT SUBINTERVAL.
+C
+      DNA=DBLE(NA)
+      DT1=DINT(T(1))
+      TEMP=DNA*T(1)
+      L=IDINT(TEMP-DT1)+1
+
+C         TC IS THE NORMALIZED CHEBYSHEV TIME (-1 .LE. TC .LE. 1)
+
+      TC=2.D0*(DMOD(TEMP,1.D0)+DT1)-1.D0
+
+C       CHECK TO SEE WHETHER CHEBYSHEV TIME HAS CHANGED,
+C       AND COMPUTE NEW POLYNOMIAL VALUES IF IT HAS.
+C       (THE ELEMENT PC(2) IS THE VALUE OF T1(TC) AND HENCE
+C       CONTAINS THE VALUE OF TC ON THE PREVIOUS CALL.)
+
+      IF(TC.NE.PC(2)) THEN
+        NP=2
+        NV=3
+        PC(2)=TC
+        TWOT=TC+TC
+      ENDIF
+C
+C       BE SURE THAT AT LEAST 'NCF' POLYNOMIALS HAVE BEEN EVALUATED
+C       AND ARE STORED IN THE ARRAY 'PC'.
+C
+      IF(NP.LT.NCF) THEN
+        DO 1 I=NP+1,NCF
+        PC(I)=TWOT*PC(I-1)-PC(I-2)
+    1   CONTINUE
+        NP=NCF
+      ENDIF
+C
+C       INTERPOLATE TO GET POSITION FOR EACH COMPONENT
+C
+      DO 2 I=1,NCM
+      PV(I,1)=0.D0
+      DO 3 J=NCF,1,-1
+      PV(I,1)=PV(I,1)+PC(J)*BUF(J,I,L)
+    3 CONTINUE
+    2 CONTINUE
+      IF(IFL.LE.1) RETURN
+C
+C       IF VELOCITY INTERPOLATION IS WANTED, BE SURE ENOUGH
+C       DERIVATIVE POLYNOMIALS HAVE BEEN GENERATED AND STORED.
+C
+      VFAC=(DNA+DNA)/T(2)
+      VC(3)=TWOT+TWOT
+      IF(NV.LT.NCF) THEN
+        DO 4 I=NV+1,NCF
+        VC(I)=TWOT*VC(I-1)+PC(I-1)+PC(I-1)-VC(I-2)
+    4   CONTINUE
+        NV=NCF
+      ENDIF
+C
+C       INTERPOLATE TO GET VELOCITY FOR EACH COMPONENT
+C
+      DO 5 I=1,NCM
+      PV(I,2)=0.D0
+      DO 6 J=NCF,2,-1
+      PV(I,2)=PV(I,2)+VC(J)*BUF(J,I,L)
+    6 CONTINUE
+      PV(I,2)=PV(I,2)*VFAC
+    5 CONTINUE
+C
+      RETURN
+C
+      END

zve/eph/pleph.f

+      SUBROUTINE PLEPH ( ET, NTARG, NCENT, RRD )
+C
+C++++++++++++++++++++++++++
+C  NOTE : Over the years, different versions of PLEPH have had a fifth argument:
+C  sometimes, an error return statement number; sometimes, a logical denoting
+C  whether or not the requested date is covered by the ephemeris.  We apologize
+C  for this inconsistency; in this present version, we use only the four necessary 
+C  arguments and do the testing outside of the subroutine.
+C
+C
+C
+C     THIS SUBROUTINE READS THE JPL PLANETARY EPHEMERIS
+C     AND GIVES THE POSITION AND VELOCITY OF THE POINT 'NTARG'
+C     WITH RESPECT TO 'NCENT'.
+C
+C     CALLING SEQUENCE PARAMETERS:
+C
+C       ET = D.P. JULIAN EPHEMERIS DATE AT WHICH INTERPOLATION
+C            IS WANTED.
+C
+C       ** NOTE THE ENTRY DPLEPH FOR A DOUBLY-DIMENSIONED TIME **
+C          THE REASON FOR THIS OPTION IS DISCUSSED IN THE 
+C          SUBROUTINE STATE
+C
+C     NTARG = INTEGER NUMBER OF 'TARGET' POINT.
+C
+C     NCENT = INTEGER NUMBER OF CENTER POINT.
+C
+C            THE NUMBERING CONVENTION FOR 'NTARG' AND 'NCENT' IS:
+C
+C                1 = MERCURY           8 = NEPTUNE
+C                2 = VENUS             9 = PLUTO
+C                3 = EARTH            10 = MOON
+C                4 = MARS             11 = SUN
+C                5 = JUPITER          12 = SOLAR-SYSTEM BARYCENTER
+C                6 = SATURN           13 = EARTH-MOON BARYCENTER
+C                7 = URANUS           14 = NUTATIONS (LONGITUDE AND OBLIQ)
+C                            15 = LIBRATIONS, IF ON EPH FILE
+C
+C             (IF NUTATIONS ARE WANTED, SET NTARG = 14. FOR LIBRATIONS,
+C              SET NTARG = 15. SET NCENT=0.)
+C
+C      RRD = OUTPUT 6-WORD D.P. ARRAY CONTAINING POSITION AND VELOCITY
+C            OF POINT 'NTARG' RELATIVE TO 'NCENT'. THE UNITS ARE AU AND
+C            AU/DAY. FOR LIBRATIONS THE UNITS ARE RADIANS AND RADIANS
+C            PER DAY. IN THE CASE OF NUTATIONS THE FIRST FOUR WORDS OF
+C            RRD WILL BE SET TO NUTATIONS AND RATES, HAVING UNITS OF
+C            RADIANS AND RADIANS/DAY.
+C
+C            The option is available to have the units in km and km/sec.
+C            For this, set km=.true. in the STCOMX common block.
+C
+
+      IMPLICIT REAL*8 (A-H,O-Z)
+
+      REAL*8 RRD(6),ET2Z(2),ET2(2),PV(6,13)
+      REAL*8 SS(3),CVAL(400),PVSUN(6),zips(2)
+C     pnut был обозначен как скаляр 08.09.2021
+      real*8 pnut(6)
+      data zips/2*0.d0/
+
+      LOGICAL BSAVE,KM,BARY
+      LOGICAL FIRST
+      DATA FIRST/.TRUE./
+
+      INTEGER LIST(12),IPT(39),DENUM
+
+
+      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,DENUM,NCON,IPT
+
+      COMMON/STCOMX/KM,BARY,PVSUN
+
+C     INITIALIZE ET2 FOR 'STATE' AND SET UP COMPONENT COUNT
+C
+      ET2(1)=ET
+      ET2(2)=0.D0
+      GO TO 11
+
+C     ENTRY POINT 'DPLEPH' FOR DOUBLY-DIMENSIONED TIME ARGUMENT 
+C          (SEE THE DISCUSSION IN THE SUBROUTINE STATE)
+
+      ENTRY DPLEPH(ET2Z,NTARG,NCENT,RRD)
+
+      ET2(1)=ET2Z(1)
+      ET2(2)=ET2Z(2)
+
+  11  IF(FIRST) CALL STATE(zips,list,pv,pnut)
+      FIRST=.FALSE.
+
+C  96  IF(NTARG .EQ. NCENT) RETURN
+      IF(NTARG .EQ. NCENT) RETURN
+
+      DO I=1,12
+      LIST(I)=0
+      ENDDO
+
+C     CHECK FOR NUTATION CALL
+
+      IF(NTARG.NE.14) GO TO 97
+        IF(IPT(35).GT.0) THEN
+          LIST(11)=2
+          CALL STATE(ET2,LIST,PV,RRD)
+          RETURN
+        ELSE
+          do i=1,4
+          rrd(i)=0.d0
+          enddo
+          WRITE(6,297)
+  297     FORMAT(' *****  NO NUTATIONS ON THE EPHEMERIS FILE  *****')
+          STOP
+        ENDIF
+
+C     CHECK FOR LIBRATIONS
+
+  97  do i=1,6
+      rrd(i)=0.d0
+      enddo
+
+      IF(NTARG.NE.15) GO TO 98
+        IF(IPT(38).GT.0) THEN
+          LIST(12)=2
+          CALL STATE(ET2,LIST,PV,RRD)
+          DO I=1,6
+          RRD(I)=PV(I,11)
+          ENDDO
+          RETURN
+        ELSE
+          WRITE(6,298)
+  298     FORMAT(' *****  NO LIBRATIONS ON THE EPHEMERIS FILE  *****')
+          STOP
+        ENDIF
+
+C       FORCE BARYCENTRIC OUTPUT BY 'STATE'
+
+  98  BSAVE=BARY
+      BARY=.TRUE.
+
+C       SET UP PROPER ENTRIES IN 'LIST' ARRAY FOR STATE CALL
+
+      DO I=1,2
+      K=NTARG
+      IF(I .EQ. 2) K=NCENT
+      IF(K .LE. 10) LIST(K)=2
+      IF(K .EQ. 10) LIST(3)=2
+      IF(K .EQ. 3) LIST(10)=2
+      IF(K .EQ. 13) LIST(3)=2
+      ENDDO
+
+C       MAKE CALL TO STATE
+
+      CALL STATE(ET2,LIST,PV,RRD)
+
+      IF(NTARG .EQ. 11 .OR. NCENT .EQ. 11) THEN
+      DO I=1,6
+      PV(I,11)=PVSUN(I)
+      ENDDO
+      ENDIF
+
+      IF(NTARG .EQ. 12 .OR. NCENT .EQ. 12) THEN
+      DO I=1,6
+      PV(I,12)=0.D0
+      ENDDO
+      ENDIF
+
+      IF(NTARG .EQ. 13 .OR. NCENT .EQ. 13) THEN
+      DO I=1,6
+      PV(I,13)=PV(I,3)
+      ENDDO
+      ENDIF
+
+      IF(NTARG*NCENT .EQ. 30 .AND. NTARG+NCENT .EQ. 13) THEN
+      DO I=1,6
+      PV(I,3)=0.D0
+      ENDDO
+      GO TO 99
+      ENDIF
+
+      IF(LIST(3) .EQ. 2) THEN
+      DO I=1,6
+      PV(I,3)=PV(I,3)-PV(I,10)/(1.D0+EMRAT)
+      ENDDO
+      ENDIF
+
+      IF(LIST(10) .EQ. 2) THEN
+      DO I=1,6
+      PV(I,10)=PV(I,3)+PV(I,10)
+      ENDDO
+      ENDIF
+
+  99  DO I=1,6
+      RRD(I)=PV(I,NTARG)-PV(I,NCENT)
+      ENDDO
+
+      BARY=BSAVE
+
+      RETURN
+      END

zve/eph/split.f

+      SUBROUTINE SPLIT(TT,FR)
+C
+C+++++++++++++++++++++++++
+C
+C     THIS SUBROUTINE BREAKS A D.P. NUMBER INTO A D.P. INTEGER
+C     AND A D.P. FRACTIONAL PART.
+C
+C     CALLING SEQUENCE PARAMETERS:
+C
+C       TT = D.P. INPUT NUMBER
+C
+C       FR = D.P. 2-WORD OUTPUT ARRAY.
+C            FR(1) CONTAINS INTEGER PART
+C            FR(2) CONTAINS FRACTIONAL PART
+C
+C            FOR NEGATIVE INPUT NUMBERS, FR(1) CONTAINS THE NEXT
+C            MORE NEGATIVE INTEGER; FR(2) CONTAINS A POSITIVE FRACTION.
+C
+C       CALLING SEQUENCE DECLARATIONS
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+
+      REAL*8 FR(2)
+
+C       MAIN ENTRY -- GET INTEGER AND FRACTIONAL PARTS
+
+      FR(1)=DINT(TT)
+      FR(2)=TT-FR(1)
+
+      IF(TT.GE.0.D0 .OR. FR(2).EQ.0.D0) RETURN
+
+C       MAKE ADJUSTMENTS FOR NEGATIVE INPUT NUMBER
+
+      FR(1)=FR(1)-1.D0
+      FR(2)=FR(2)+1.D0
+
+      RETURN
+
+      END

zve/eph/state.f

+      SUBROUTINE STATE(ET2,LIST,PV,PNUT)
+C
+C++++++++++++++++++++++++++++++++
+C
+C THIS SUBROUTINE READS AND INTERPOLATES THE JPL PLANETARY EPHEMERIS FILE
+C
+C     CALLING SEQUENCE PARAMETERS:
+C
+C     INPUT:
+C
+C         ET2   DP 2-WORD JULIAN EPHEMERIS EPOCH AT WHICH INTERPOLATION
+C               IS WANTED.  ANY COMBINATION OF ET2(1)+ET2(2) WHICH FALLS
+C               WITHIN THE TIME SPAN ON THE FILE IS A PERMISSIBLE EPOCH.
+C
+C                A. FOR EASE IN PROGRAMMING, THE USER MAY PUT THE
+C                   ENTIRE EPOCH IN ET2(1) AND SET ET2(2)=0.
+C
+C                B. FOR MAXIMUM INTERPOLATION ACCURACY, SET ET2(1) =
+C                   THE MOST RECENT MIDNIGHT AT OR BEFORE INTERPOLATION
+C                   EPOCH AND SET ET2(2) = FRACTIONAL PART OF A DAY
+C                   ELAPSED BETWEEN ET2(1) AND EPOCH.
+C
+C                C. AS AN ALTERNATIVE, IT MAY PROVE CONVENIENT TO SET
+C                   ET2(1) = SOME FIXED EPOCH, SUCH AS START OF INTEGRATION,
+C                   AND ET2(2) = ELAPSED INTERVAL BETWEEN THEN AND EPOCH.
+C
+C        LIST   12-WORD INTEGER ARRAY SPECIFYING WHAT INTERPOLATION
+C               IS WANTED FOR EACH OF THE BODIES ON THE FILE.
+C
+C                         LIST(I)=0, NO INTERPOLATION FOR BODY I
+C                                =1, POSITION ONLY
+C                                =2, POSITION AND VELOCITY
+C
+C               THE DESIGNATION OF THE ASTRONOMICAL BODIES BY I IS:
+C
+C                         I = 1: MERCURY
+C                           = 2: VENUS
+C                           = 3: EARTH-MOON BARYCENTER
+C                           = 4: MARS
+C                           = 5: JUPITER
+C                           = 6: SATURN
+C                           = 7: URANUS
+C                           = 8: NEPTUNE
+C                           = 9: PLUTO
+C                           =10: GEOCENTRIC MOON
+C                           =11: NUTATIONS IN LONGITUDE AND OBLIQUITY
+C                           =12: LUNAR LIBRATIONS (IF ON FILE)
+C
+C
+C     OUTPUT:
+C
+C          PV   DP 6 X 11 ARRAY THAT WILL CONTAIN REQUESTED INTERPOLATED
+C               QUANTITIES.  THE BODY SPECIFIED BY LIST(I) WILL HAVE ITS
+C               STATE IN THE ARRAY STARTING AT PV(1,I).  (ON ANY GIVEN
+C               CALL, ONLY THOSE WORDS IN 'PV' WHICH ARE AFFECTED BY THE
+C               FIRST 10 'LIST' ENTRIES (AND BY LIST(12) IF LIBRATIONS ARE
+C               ON THE FILE) ARE SET.  THE REST OF THE 'PV' ARRAY
+C               IS UNTOUCHED.)  THE ORDER OF COMPONENTS STARTING IN
+C               PV(1,I) IS: X,Y,Z,DX,DY,DZ.
+C
+C               ALL OUTPUT VECTORS ARE REFERENCED TO THE EARTH MEAN
+C               EQUATOR AND EQUINOX OF J2000 IF THE DE NUMBER IS 200 OR
+C               GREATER; OF B1950 IF THE DE NUMBER IS LESS THAN 200. 
+C
+C               THE MOON STATE IS ALWAYS GEOCENTRIC; THE OTHER NINE STATES 
+C               ARE EITHER HELIOCENTRIC OR SOLAR-SYSTEM BARYCENTRIC, 
+C               DEPENDING ON THE SETTING OF COMMON FLAGS (SEE BELOW).
+C
+C               LUNAR LIBRATIONS, IF ON FILE, ARE PUT INTO PV(K,11) IF
+C               LIST(12) IS 1 OR 2.
+C
+C         NUT   DP 4-WORD ARRAY THAT WILL CONTAIN NUTATIONS AND RATES,
+C               DEPENDING ON THE SETTING OF LIST(11).  THE ORDER OF
+C               QUANTITIES IN NUT IS:
+C
+C                        D PSI  (NUTATION IN LONGITUDE)
+C                        D EPSILON (NUTATION IN OBLIQUITY)
+C                        D PSI DOT
+C                        D EPSILON DOT
+C
+C           *   STATEMENT # FOR ERROR RETURN, IN CASE OF EPOCH OUT OF
+C               RANGE OR I/O ERRORS.
+C
+C
+C     COMMON AREA STCOMX:
+C
+C          KM   LOGICAL FLAG DEFINING PHYSICAL UNITS OF THE OUTPUT
+C               STATES. KM = .TRUE., KM AND KM/SEC
+C                          = .FALSE., AU AND AU/DAY
+C               DEFAULT VALUE = .FALSE.  (KM DETERMINES TIME UNIT
+C               FOR NUTATIONS AND LIBRATIONS.  ANGLE UNIT IS ALWAYS RADIANS.)
+C
+C        BARY   LOGICAL FLAG DEFINING OUTPUT CENTER.
+C               ONLY THE 9 PLANETS ARE AFFECTED.
+C                        BARY = .TRUE. =\ CENTER IS SOLAR-SYSTEM BARYCENTER
+C                             = .FALSE. =\ CENTER IS SUN
+C               DEFAULT VALUE = .FALSE.
+C
+C       PVSUN   DP 6-WORD ARRAY CONTAINING THE BARYCENTRIC POSITION AND
+C               VELOCITY OF THE SUN.
+C
+C
+      IMPLICIT REAL*8 (A-H,O-Z)
+
+      SAVE
+C     Было PV(6,12) иÑправлено PV(6,13) 08.09.2021
+C     Было PNUT(4) иÑправлено PNUT(6) 08.09.2021
+      REAL*8 ET2(2),PV(6,13),PNUT(6),T(2),PJD(4),BUF(1500),
+     . SS(3),CVAL(400),PVSUN(6)
+
+      INTEGER LIST(12),IPT(3,13)
+
+      LOGICAL FIRST
+      DATA FIRST/.TRUE./
+
+      CHARACTER*6 TTL(14,3),CNAM(400)
+      CHARACTER*80 NAMFIL
+
+      LOGICAL KM,BARY
+*  â êèëîìåòðàõ:
+      DATA KM/.TRUE./
+
+      COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
+      COMMON/CHRHDR/CNAM,TTL
+      COMMON/STCOMX/KM,BARY,PVSUN
+
+C
+C       ENTRY POINT - 1ST TIME IN, GET POINTER DATA, ETC., FROM EPH FILE
+C
+      IF(FIRST) THEN
+        FIRST=.FALSE.
+
+C ************************************************************************
+C ************************************************************************
+
+C THE USER MUST SELECT ONE OF THE FOLLOWING BY DELETING THE 'C' IN COLUMN 1
+
+C ************************************************************************
+
+C        CALL FSIZER1(NRECL,KSIZE,NRFILE,NAMFIL)
+C        CALL FSIZER2(NRECL,KSIZE,NRFILE,NAMFIL)
+        CALL FSIZER3(NRECL,KSIZE,NRFILE,NAMFIL)
+
+      IF(NRECL .EQ. 0) WRITE(*,*)'  ***** FSIZER IS NOT WORKING *****'
+
+C ************************************************************************
+C ************************************************************************
+
+      IRECSZ=NRECL*KSIZE
+      NCOEFFS=KSIZE/2
+
+        OPEN(NRFILE,
+     *       FILE=NAMFIL,
+     *       ACCESS='DIRECT',
+     *       FORM='UNFORMATTED',
+     *       RECL=IRECSZ,
+     *       STATUS='OLD')
+
+      READ(NRFILE,REC=1)TTL,CNAM,SS,NCON,AU,EMRAT,
+     . ((IPT(I,J),I=1,3),J=1,12),NUMDE,(IPT(I,13),I=1,3)
+
+      READ(NRFILE,REC=2)CVAL
+
+      NRL=0
+
+      ENDIF
+
+
+C       ********** MAIN ENTRY POINT **********
+
+
+      IF(ET2(1) .EQ. 0.D0) RETURN
+
+      S=ET2(1)-.5D0
+      CALL SPLIT(S,PJD(1))
+      CALL SPLIT(ET2(2),PJD(3))
+      PJD(1)=PJD(1)+PJD(3)+.5D0
+      PJD(2)=PJD(2)+PJD(4)
+      CALL SPLIT(PJD(2),PJD(3))
+      PJD(1)=PJD(1)+PJD(3)
+
+C       ERROR RETURN FOR EPOCH OUT OF RANGE
+
+      IF(PJD(1)+PJD(4).LT.SS(1) .OR. PJD(1)+PJD(4).GT.SS(2)) GO TO 98
+
+C       CALCULATE RECORD # AND RELATIVE TIME IN INTERVAL
+
+      NR=IDINT((PJD(1)-SS(1))/SS(3))+3
+      IF(PJD(1).EQ.SS(2)) NR=NR-1
+
+        tmp1 = DBLE(NR-3)*SS(3) + SS(1)
+        tmp2 = PJD(1) - tmp1
+        T(1) = (tmp2 + PJD(4))/SS(3)
+
+C       READ CORRECT RECORD IF NOT IN CORE
+
+      IF(NR.NE.NRL) THEN
+        NRL=NR
+        READ(NRFILE,REC=NR,ERR=99)(BUF(K),K=1,NCOEFFS)
+      ENDIF
+
+      IF(KM) THEN
+      T(2)=SS(3)*86400.D0
+      AUFAC=1.D0
+      ELSE
+      T(2)=SS(3)
+      AUFAC=1.D0/AU
+      ENDIF
+
+C   INTERPOLATE SSBARY SUN
+
+      CALL INTERP(BUF(IPT(1,11)),T,IPT(2,11),3,IPT(3,11),2,PVSUN)
+
+      DO I=1,6
+      PVSUN(I)=PVSUN(I)*AUFAC
+      ENDDO
+
+C   CHECK AND INTERPOLATE WHICHEVER BODIES ARE REQUESTED
+
+      DO 4 I=1,10
+      IF(LIST(I).EQ.0) GO TO 4
+
+      CALL INTERP(BUF(IPT(1,I)),T,IPT(2,I),3,IPT(3,I),
+     & LIST(I),PV(1,I))
+
+      DO J=1,6
+       IF(I.LE.9 .AND. .NOT.BARY) THEN
+       PV(J,I)=PV(J,I)*AUFAC-PVSUN(J)
+       ELSE
+       PV(J,I)=PV(J,I)*AUFAC
+       ENDIF
+      ENDDO
+
+   4  CONTINUE
+
+C       DO NUTATIONS IF REQUESTED (AND IF ON FILE)
+
+      IF(LIST(11).GT.0 .AND. IPT(2,12).GT.0)
+     * CALL INTERP(BUF(IPT(1,12)),T,IPT(2,12),2,IPT(3,12),
+     * LIST(11),PNUT)
+
+C       GET LIBRATIONS IF REQUESTED (AND IF ON FILE)
+
+      IF(LIST(12).GT.0 .AND. IPT(2,13).GT.0)
+     * CALL INTERP(BUF(IPT(1,13)),T,IPT(2,13),3,IPT(3,13),
+     * LIST(12),PV(1,11))
+
+      RETURN
+
+  98  WRITE(0,198)ET2(1)+ET2(2),SS(1),SS(2)
+ 198  format(' ***  Requested JED,',f12.2,
+     * ' not within ephemeris limits,',2f12.2,'  ***')
+
+      stop
+
+   99 WRITE(*,'(2F12.2,A80)')ET2,'ERROR RETURN IN STATE'
+
+      STOP
+
+      END

zve/gencatalog/.gitignore

+write_buildinfo.f95
+!Doxyfile
\ No newline at end of file

zve/gencatalog/makefile

+MAIN=gencatalog
+include ../makefile.inc
\ No newline at end of file

zve/makefile.inc

+PLATFORM:=$(shell uname -s |awk -F- '{print $$1}')
+
+FC=gfortran
+RM=rm
+
+BIN = ../../bin
+BUILDDIR = ../../build/$(PLATFORM)
+MODDIR = $(BUILDDIR)/mod
+LIBDIR = $(BUILDDIR)/lib
+SHAREDIR = $(BUILDDIR)/share
+DOCDIR = $(BUILDDIR)/docs/$(MAIN)
+LATEXDIR = $(DOCDIR)/latex
+
+LDFLAGS  = -lrdjson -lfson -lsolv -lsofa -lmemc_gate -lgravity -lmemcached -lsofa_c -lstdc++ -L$(LIBDIR)
+FCFLAGS  = -MD -Wall -fmax-stack-var-size=3932160 -ffree-line-length-none -J$(MODDIR) -cpp
+
+ifeq ($(PLATFORM),MINGW32_NT)
+    EXE=.exe
+else
+    EXE=
+    FCFLAGS += -fPIC
+endif
+
+ifeq ($(PLATFORM), FreeBSD)
+    RPATH_FREEBSD=$(shell gcc -print-search-dirs | grep libraries | tr -d "= " | awk -F: '{print $2}' | awk -F'/gcc/' '{print $1}')
+    LDFLAGS += -lexecinfo -rpath $(RPATH_FREEBSD)
+endif
+
+SRC = $(wildcard *.f*)
+OBJS = $(patsubst %.f,%.o,$(patsubst %.f95,%.o,$(wildcard *.f*)))
+DEPS = $(patsubst %.f,%.d,$(patsubst %.f95,%.d,$(wildcard *.f*)))
+
+LIBDEPS = libsofa.a libsolv.a libmemc_gate.a libgravity.a librdjson.a libfson.a
+LIBS = $(addprefix $(LIBDIR)/, $(LIBDEPS))
+
+RELDIR = $(BUILDDIR)/obj/solvers/$(MAIN)
+RELEXE = $(BUILDDIR)/bin/$(MAIN)$(EXE)
+RELOBJS = $(addprefix $(RELDIR)/, $(OBJS))
+RELDEPS = $(addprefix $(RELDIR)/, $(DEPS))
+
+ifdef DEBUG
+	RELCFLAGS = -g -O0 -DDEBUG=$(DEBUG)
+else
+	RELCFLAGS = -O2
+endif
+
+LAST = $(shell ls -atR *.f* | head -n 1)
+COMMIT = $(shell git rev-parse HEAD)
+
+all: dirs release
+
+.PHONY: dirs
+dirs:
+	mkdir -p $(MODDIR) $(SHAREDIR) $(LIBDIR) $(DEPSDIR)
+	mkdir -p $(RELDIR) $(BUILDDIR)/bin
+
+.PHONY: buildinfo
+buildinfo:
+	if [ "$(LAST)" != "write_buildinfo.f95" ]; then ../buildinfo; fi
+
+write_buildinfo.f95: $(LIBS)
+	../buildinfo
+	$(FC) -c $(FCFLAGS) $(RELCFLAGS) -o $(RELDIR)/write_buildinfo.o write_buildinfo.f95
+
+release: buildinfo $(RELEXE)
+$(RELEXE): $(RELOBJS) $(LIBS)
+	$(FC) -o $(RELEXE) $^ $(LDFLAGS)
+$(RELDIR)/%.o: %.f*
+	$(FC) -c $(FCFLAGS) $(RELCFLAGS) -o $@ $<
+include $(wildcard $(RELDIR)/*.d)
+
+include $(BIN)/docs.mk
+
+.PHONY: transform
+
+transform:
+	$(BIN)/each
+
+.PHONY: clean
+
+clean:
+	$(RM) -rf $(RELOBJS) $(RELEXE) $(DBGOBJS) $(DBGEXE) $(RELDEPS) $(DBGDEPS) $(DOCDIR)
+
+
+

zve/math/.gitignore

+!Doxyfile
\ No newline at end of file

zve/math/Linterpolation.f95

+!-----------------------------------------------------------------------
+!>  @author Нью-Ком Технолоджис
+!>  @file Linterpolation.f95
+!>  @date Mar 6, 2019
+!>  @addtogroup math
+!>  @{
+!-----------------------------------------------------------------------
+      SUBROUTINE Linterpolation (T,Teph,Yeph,Y)
+!-----------------------------------------------------------------------
+!
+!> @brief    Подпрограмма интерполяции методом Лагранжа 
+!
+!> @details  Получение интерполяционного значения методом Лагранжа по 8 точкам 
+!
+!-----------------------------------------------------------------------
+!>  Входные параметры
+!> @param[in]  T               R*8 - Значение аргумента  
+!> @param[in]  Teph(8)         R*8 - Массив значений аргумента в узлах интерполяции
+!> @param[in]  Yeph(8)         R*8 - Массив значений функций в узлах интерполяции
+!>  Выходные параметры
+!> @param[out] Y               R*8 - Искомое значение функции.
+!-----------------------------------------------------------------------
+
+    use const_m
+    use adaps_m
+
+    IMPLICIT REAL*8 (A-H,O-Z)
+
+    REAL(8) :: T          ! Значение аргумента
+    REAL(8) :: Teph(8)    ! Массив значений аргумента в узлах интерполяции
+    REAL(8) :: Yeph(8)    ! Массив значений функций в узлах интерполяции
+    REAL(8) :: Y          ! искомое значение функции
+    INTEGER :: i,j        ! переменные цикла
+    
+    REAL(8) :: L          ! базисный полином
+
+    Y=0.0D00              ! Инициализация
+    do i=1,8
+        ! определяем базисный полином
+        L=1.0
+        do j=1,8
+            if(i.ne.j) L=L*(T-Teph(j))/(Teph(i)-Teph(j))
+        enddo
+        Y=Y+Yeph(i)*L
+    enddo    
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "T",T
+    write(ERROR_UNIT,*) "Teph",Teph
+    write(ERROR_UNIT,*) "Yeph",Yeph
+    write(ERROR_UNIT,*) "rezult",Y
+#endif    
+
+    RETURN
+END

zve/math/diffANGL.f95

+!-----------------------------------------------------------------------
+!
+!  Подпрограмма определения разности двух углов с учетом четвертей
+!  diffANGL = ALF - ALF0
+!   (вспомогательная)
+!
+!
+!-----------------------------------------------------------------------
+!
+!  Входные параметры
+!
+!    ALF             R*8 - угол измеренный
+!    ALF0            R*8 - угол прогнозный
+!-----------------------------------------------------------------------
+
+REAL*8 FUNCTION diffANGL(ALF,ALF0)  
+
+    IMPLICIT NONE
+
+    REAL(8) :: ALF                  ! угол
+    real(8) :: ALF0                 ! угол
+    real(8) :: D2PI = 6.283185307179586476925287D0
+    REAL(8) :: PI4 =  0.785398163397448309615661D0
+
+    if(ALF0.GT.D2PI-PI4.AND.ALF.LT.PI4) then         
+        diffANGL = (ALF+D2PI)-ALF0
+    ELSEIF(ALF.GT.D2PI-PI4.AND.ALF0.LT.PI4) then
+        diffANGL=ALF-(ALF0+D2PI)
+    ELSE
+        diffANGL=ALF-ALF0 
+    ENDIF
+
+    return
+end

zve/math/dpolrt.f

+      SUBROUTINE DPOLRT(XCOF,COF,M,ROOTR,ROOTI,IER)
+      DIMENSION XCOF(1),COF(1),ROOTR(1),ROOTI(1)
+      DOUBLE PRECISION XCOF,COF,ROOTR,ROOTI,XO,YO,X,Y,XPR,YPR,UX,UY,V,
+     1YT,XT,U,XT2,YT2,SUMSQ,DX,DY,TEMP,ALPHA
+      IFIT=0
+      YPR=0.0D00
+      XPR=0.0D00
+      N=M
+      IER=0
+      IF(XCOF(N+1))10,25,10
+   10 IF(N) 15,15,35
+   15 IER=1
+   20 RETURN
+   25 IER=2
+      GO TO 20
+   35 NX=N
+      NXX=N+1
+      N2=1
+      KJ1=N+1
+      DO 40 L=1,KJ1
+      MT=KJ1-L+1
+   40 COF(MT)=XCOF(L)
+   45 XO=.500101D-2
+      YO=.1000101D-1
+      IN=0
+   50 X=XO
+      XO=-10.D0 *YO
+      YO=-10.D0 *X
+      X=XO
+      Y=YO
+      IN=IN+1
+      GO TO 59
+   55 IFIT=1
+      XPR=X
+      YPR=Y
+   59 ICT=0
+   60 UX=0.D0
+      UY=0.D0
+      V=0.D0
+      YT=0.D0
+      XT=1.D0
+      U=COF(N+1)
+      IF(U) 65,130,65
+   65 DO 70 I=1,N
+      L=N-I+1
+      TEMP=COF(L)
+      XT2=X*XT-Y*YT
+      YT2=X*YT+Y*XT
+      U=U+TEMP*XT2
+      V=V+TEMP*YT2
+      FI=I
+      UX=UX+FI*XT*TEMP
+      UY=UY-FI*YT*TEMP
+      XT=XT2
+   70 YT=YT2
+      SUMSQ=UX*UX+UY*UY
+      IF(SUMSQ) 75,110,75
+   75 DX=(V*UY-U*UX)/SUMSQ
+      X=X+DX
+      DY=-(U*UY+V*UX)/SUMSQ
+      Y=Y+DY
+      IF(DABS(DY)+DABS(DX)-1.0D-12) 100,80,80
+   80 ICT=ICT+1
+      IF(ICT-500) 60,85,85
+   85 IF(IFIT)100,90,100
+   90 IF(IN-5) 50,95,95
+   95 IER=3
+      GO TO 20
+  100 DO 105 L=1,NXX
+      MT=KJ1-L+1
+      TEMP=XCOF(MT)
+      XCOF(MT)=COF(L)
+  105 COF(L)=TEMP
+      ITEMP=N
+      N=NX
+      NX=ITEMP
+      IF(IFIT) 120,55,120
+  110 IF(IFIT) 115,50,115
+  115 X=XPR
+      Y=YPR
+  120 IFIT=0
+      IF(DABS(Y)-1.0D-10*DABS(X)) 135,125,125
+  125 ALPHA=X+X
+      SUMSQ=X*X+Y*Y
+      N=N-2
+      GO TO 140
+  130 X=0.D0
+      NX=NX-1
+      NXX=NXX-1
+  135 Y=0.D0
+      SUMSQ=0.D0
+      ALPHA=X
+      N=N-1
+  140 COF(2)=COF(2)+ALPHA*COF(1)
+      DO 150 L=2,N
+  150 COF(L+1)=COF(L+1)+ALPHA*COF(L)-SUMSQ*COF(L-1)
+  155 ROOTI(N2)=Y
+      ROOTR(N2)=X
+      N2=N2+1
+      IF(SUMSQ)160,165,160
+  160 Y=-Y
+      SUMSQ=0.D0
+      GO TO 155
+  165 IF(N) 20,20,45
+      END

zve/math/nc_dminv.f95

+!>  @file nc_dminv.f95
+!>  @date Jun 18, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы
+!> @details  Подпрограмма выполняет обращение матрицы
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[inout]   A(N,N)  R*8 - обращаемая матрица A (туда же списывается результат обращения)
+!> @param[in]      N       I*4 - размерность матрицы А
+!> @param[out]     D       R*4 - дискриминант
+!> @param[inout]   L(N)    I*4 - вспомогательный целочисленный массив
+!> @param[inout]   M(N)    I*4 - вспомогательный целочисленный массив
+!>
+SUBROUTINE nc_dminv(A,N,D,L,M)
+
+    implicit none
+
+    real(8), intent(inout) :: A(1)  !< обращаемая матрица
+    integer, intent(in)    :: N     !< размерность матрицы
+    real(8), intent(out)   :: D     !< дискриминатнт
+    integer, intent(inout) :: L(1)  !< вспомогательный целочисленный массив
+    integer, intent(inout) :: M(1)  !< вспомогательный целочисленный массив
+    
+    real(8) :: BIGA
+    real(8) :: HOLD
+    
+    integer :: i,j,k
+    integer :: NK,IZ,IJ,IK
+    integer :: JR,JQ,JP,JK,JI
+    integer :: KJ,KK,KI
+
+    D=1.0D00
+    NK=-N
+    
+    do k=1,N
+        NK=NK+N
+        L(K)=K
+        M(K)=K
+        KK=NK+K
+        BIGA=A(KK)
+
+        do j=k,N
+            IZ=N*(j-1)
+            do i=K,N
+                IJ=IZ+I
+                if((DABS(BIGA)-DABS(A(IJ))).lt.0.0D00) then
+                    BIGA=A(IJ)
+                    L(k)=I
+                    M(k)=J
+                endif
+            enddo
+        enddo
+
+        j=L(K)
+
+        if((j-K).gt.0.0D00) then
+            KI=K-N
+            do i=1,N
+                KI=KI+N
+                HOLD=-A(KI)
+                JI=KI-k+j
+                A(KI)=A(JI)
+                A(JI)=HOLD
+            enddo
+        endif
+
+        i=M(K)
+ 
+        IF((i-k).gt.0.0D00) then
+            JP=N*(i-1)
+            do j=1,N
+                JK=NK+j
+                JI=JP+j
+                HOLD=-A(JK)
+                A(JK)=A(JI)
+                A(JI)=HOLD
+            enddo
+        endif
+
+        if(BIGA.eq.0) then
+            D=0.D0
+            RETURN
+        endif
+
+        DO i=1,N
+            IF((I-K).ne.0) then
+                IK=NK+I
+                A(IK)=A(IK)/(-BIGA)
+            endif
+        enddo
+
+        do i=1,N
+            IK=NK+i
+            HOLD=A(IK)
+            IJ=I-N
+            do j=1,N
+                IJ=IJ+N
+                if((I-K).ne.0) then
+                    if((J-K).ne.0) then
+                        KJ=IJ-I+K
+                        A(IJ)=HOLD*A(KJ)+A(IJ)
+                    endif
+                endif
+            enddo
+        enddo    
+
+        KJ=k-N
+
+        do j=1,N
+            KJ=KJ+N
+            IF((J-K).ne.0) then
+                A(KJ)=A(KJ)/BIGA
+            endif
+        enddo
+
+        D=D*BIGA
+        A(KK)=1.D0/BIGA
+    enddo
+
+    K=N
+
+    do k=N-1,1,-1
+        i=L(k)
+        if((i-k).gt.0) then
+            JQ=N*(k-1)
+            JR=N*(i-1)
+
+            do j=1,N
+                JK=JQ+j
+                HOLD=A(JK)
+                JI=JR+j
+                A(JK)=-A(JI)
+                A(JI)=HOLD
+            enddo
+        endif
+
+        j=M(k)
+
+        if((j-k).gt.0) then
+            KI=k-n
+            do i=1,N
+                KI=KI+N
+                HOLD=A(KI)
+                JI=KI-k+j
+                A(KI)=-A(JI)
+                A(JI)=HOLD
+            enddo
+        endif
+    enddo
+
+    return
+end

zve/math/nc_dover.f

+*-----------------------------------------------------------------------
+*>  @author Нью-Ком Технолоджис 
+*>  @file nc_dover.f
+*>  @date Apr 3, 2019
+*>  @addtogroup library
+*>  @{
+*-----------------------------------------------------------------------
+      SUBROUTINE nc_dover ( PVER,M,G1,G2 )
+*-----------------------------------------------------------------------
+*>  @brief  Программа расчета верхней и нижней границы доверительного интервала
+*
+*  Входные параметры
+* @param[in]   PVER             R*8 - заданная вероятность
+*                                    (времено фиксированная 0.99)
+* @param[in]   M                I   - степень свободы (N-n-1)
+*  Выходные параметры
+* @param[out]  G1               R*8 - нижняя граница
+* @param[out]  G2               R*8 - верхняя граница
+*----------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      integer M                              ! степень свободы
+      REAL*8  PVER                           ! заданная вероятность
+      REAL*8  G1                             ! нижняя граница
+      REAL*8  G2                             ! верхняя граница
+      integer i                              ! рабочая
+
+      REAL*8 TAB(37,2)/                      ! заданные табличные значени
+     #  0.356,0.434,0.483,0.519,0.546,0.569,0.588,0.604,0.618,0.630,
+     #  0.641,0.651,0.660,0.669,0.676,0.683,0.690,0.696,0.792,0.707,
+     #  0.712,0.717,0.722,0.726,0.730,0.734,0.737,0.741,0.744,0.748,
+     #  0.774,0.793,0.808,0.820,0.829,0.838,0.845,
+     #  159.0,14.10,6.470,4.390,3.480,2.980,2.660,2.440,2.277,2.154,
+     #  2.056,1.976,1.910,1.854,1.806,1.764,1.727,1.695,1.666,1.640,
+     #  1.617,1.595,1.576,1.558,1.541,1.526,1.512,1.499,1.487,1.475,
+     #  1.390,1.336,1.299,1.272,1.259,1.233,1.219/
+
+      G1=0.0D00
+      G2=0.0D00
+      if(M.le.0) then
+         write(0,'("M=0",3I10)') PVER,M
+         return
+      endif
+      IF(M.ge.1.and.M.le.30)  then
+        G1=TAB(M,1)
+        G2=TAB(M,2)
+      else if (M.ge.31.and.M.le.100) then
+        i=int(M/10.0D00)
+        i=30+i-3
+        G1=TAB(i,1)
+        G2=TAB(i,2)
+      else
+        G1=TAB(37,1)
+        G2=TAB(37,2)
+      endif
+      return
+      end

zve/math/nc_mabc.f95

+!>  @file nc_mabc.f95
+!>  @date Jun 21, 2020
+!>
+!> @brief  Подпрограмма перемножения матриц
+!> @details  Подпрограмма выполняет перемножение двух матриц 
+!> R(NxL)=A(NxM)*B(MxL)
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   N       I*4 - количество строк в матрице А
+!> @param[in]   M       I*4 - количество строк в матрице B и столбцов в матрице A
+!> @param[in]   L       I*4 - количество столбцов в матрице B
+!> @param[in]  A(N,M)   R*8 - матрица A
+!> @param[in]  B(M,L)   R*8 - матрица B
+!> @param[out] R(N,L)   R*8 - произведение матриц A и B
+!>
+
+SUBROUTINE nc_mabc(N,M,L,A,B,R)
+
+    use adaps_m
+    
+    IMPLICIT NONE
+    
+    integer, intent(in)  :: N       !< количество строк в матрице А
+    integer, intent(in)  :: M       !< количество строк в матрице B и столбцов в матрице A
+    integer, intent(in)  :: L       !< количество столбцов в матрице B
+    REAL(8), intent(in)  :: A(N,M)  !< матрица A
+    REAL(8), intent(in)  :: B(M,L)  !< матрица B
+    REAL(8), intent(out) :: R(N,L)  !< матрица R = A*B
+
+    REAL(8) :: PR      ! рабочая переменная 
+    integer :: i,j,k   ! переменные цикла
+
+    do i=1,N
+        do k=1,L
+            PR=0.0D00
+            do j=1,M
+                PR=PR+A(i,j)*B(j,k)
+            enddo
+            R(i,k)=PR
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_mabtc.f95

+!>  @file nc_mabtc.f95
+!>  @date Jun 22, 2020
+!>
+!> @brief  Подпрограмма перемножения матриц
+!> @details  Подпрограмма выполняет перемножение двух матриц 
+!> R(NxL)=A(NxM)*BT(MxL). Матрица B транспонируется
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   N       I*4 - количество строк в матрице А
+!> @param[in]   M       I*4 - количество строк в матрице B и столбцов в матрице A
+!> @param[in]   L       I*4 - количество столбцов в матрице B
+!> @param[in]  A(N,M)   R*8 - матрица A
+!> @param[in]  B(L,M)   R*8 - матрица B
+!> @param[out] R(N,L)   R*8 - произведение матриц A и B
+!>
+SUBROUTINE nc_mabtc(N,M,L,A,B,R)
+
+    use adaps_m
+    
+    IMPLICIT NONE
+    
+    integer, intent(in)  :: N       !< количество строк в матрице А
+    integer, intent(in)  :: M       !< количество строк в матрице B и столбцов в матрице A
+    integer, intent(in)  :: L       !< количество столбцов в матрице B
+    REAL(8), intent(in)  :: A(N,M)  !< матрица A
+    REAL(8), intent(in)  :: B(L,M)  !< матрица B
+    REAL(8), intent(out) :: R(N,L)  !< матрица R = A*B
+
+    REAL(8) :: PR      ! рабочая переменная 
+    integer :: i,j,k   ! переменные цикла
+
+    do I=1,N
+        do K=1,L
+            PR=0.0D00
+            do J=1,M
+                PR=PR+A(I,J)*B(K,J)
+            enddo
+            R(I,K)=PR
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_matbc.f95

+!>  @file nc_matbc.f95
+!>  @date Jun 22, 2020
+!>
+!> @brief  Подпрограмма перемножения матриц
+!> @details  Подпрограмма выполняет перемножение двух матриц 
+!> R(NxL)=AT(NxM)*B(MxL). Матрица A транспонируется
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   N       I*4 - количество строк в матрице А
+!> @param[in]   M       I*4 - количество строк в матрице B и столбцов в матрице A
+!> @param[in]   L       I*4 - количество столбцов в матрице B
+!> @param[in]  A(M,N)   R*8 - матрица A
+!> @param[in]  B(M,L)   R*8 - матрица B
+!> @param[out] R(N,L)   R*8 - произведение матриц A и B
+!>
+SUBROUTINE nc_matbc(N,M,L,A,B,R)
+
+    use adaps_m
+    
+    IMPLICIT NONE
+    
+    integer, intent(in)  :: N       !< количество строк в матрице А
+    integer, intent(in)  :: M       !< количество строк в матрице B и столбцов в матрице A
+    integer, intent(in)  :: L       !< количество столбцов в матрице B
+    REAL(8), intent(in)  :: A(M,N)  !< матрица A
+    REAL(8), intent(in)  :: B(M,L)  !< матрица B
+    REAL(8), intent(out) :: R(N,L)  !< матрица R = A*B
+
+    REAL(8) :: PR      ! рабочая переменная 
+    integer :: i,j,k   ! переменные цикла
+
+    DO I=1,N
+        DO K=1,L
+            PR=0.0D00
+            DO J=1,M
+                PR=PR+A(J,I)*B(J,K)
+            enddo
+            R(I,K)=PR
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_mdet2.f95

+!>  @file nc_mdet2.f95
+!>  @date Jun 28, 2020
+!>
+!> @brief  Расчет определителя матрицы  2x2
+!> @details  Подпрограмма выполняет расчет определителя матрицы 2x2
+!>
+!> Входные пapaмeтpы: 
+!>
+!> @param[in]    A(2,2)    R*8 - матрица A
+!>
+!> Возвращаемое значение:
+!>
+!>               nc_mdet2  R*8 - определитель матрицы
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+real(8) FUNCTION nc_mdet2 ( A )
+
+    IMPLICIT NONE
+
+    real(8), intent(in) :: A(2,2)
+
+    nc_mdet2 = A(1,1)*A(2,2)-A(1,2)*A(2,1)
+
+    return
+end

zve/math/nc_mdet3.f95

+!>  @file nc_mdet3.f95
+!>  @date Jun 26, 2020
+!>
+!> @brief  Расчет определителя матрицы  3x3
+!> @details  Подпрограмма выполняет расчет определителя матрицы 3x3
+!>
+!> Входные пapaмeтpы: 
+!>
+!> @param[in]    A(3,3)    R*8 - матрица A
+!>
+!> Возвращаемое значение:
+!>
+!>               nc_mdet3  R*8 - определитель матрицы
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+real(8) FUNCTION nc_mdet3 ( A )
+
+    IMPLICIT NONE
+
+    real(8), intent(in) :: A(3,3)
+
+    nc_mdet3 = A(1,1)*A(2,2)*A(3,3)+&
+               A(1,2)*A(2,3)*A(3,1)+&
+               A(1,3)*A(2,1)*A(3,2)-&
+               A(1,3)*A(2,2)*A(3,1)-&
+               A(1,2)*A(2,1)*A(3,3)-&
+               A(1,1)*A(2,3)*A(3,2)
+
+     return
+end

zve/math/nc_minv2.f95

+!>  @file nc_minv2.f95
+!>  @date Jun 28, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы 2x2
+!> @details  Подпрограмма выполняет обращение матрицы 2x2
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]    MA(2,2)   R*8 - обращаемая матрица A
+!> @param[out]   AOBR(2,2) R*8 - искомая обратная матрица
+!> @param[out]   RC        I*4 - код возврата
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+SUBROUTINE nc_minv2 ( MA, AOBR, RC )
+
+    use adaps_m
+
+    implicit none
+
+    real(8), intent(in)   :: MA(2,2)   !< обращаемая матрица
+    real(8), intent(out)  :: AOBR(2,2) !< обратная матрица
+    integer, intent(out)  :: RC        !< код возврата
+
+    real(8) :: nc_mdet2     ! функция поиска определителя матрицы 3x3
+
+    real(8) :: A(2,2)       ! Нормированная матрица
+    real(8) :: DET          ! определитель матрицы 3x3
+    real(8) :: S            ! множитель
+    
+    integer :: i,j          ! переменные цикла
+
+    RC = 0
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv2 : MA:"
+    write(ERROR_UNIT,'(2G20.12)') ((MA(i,j),j=1,2),i=1,2)
+#endif
+
+    S=0.0D00
+    do i=1,2
+        do j=1,2
+            IF(DABS(MA(I,J)).gt.0.1D-50) S=S+DLOG(DABS(MA(i,j)))
+        enddo
+    enddo
+
+    S=S/4
+    S=1.0D00/DEXP(S)
+    
+    A(:,:)=MA(:,:)*S
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv2 : A:"
+    write(ERROR_UNIT,'(2G20.12)') ((A(i,j),j=1,2),i=1,2)
+#endif
+
+    ! Находим определитель 
+    DET = nc_mdet2(A)
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv2 : DET=",DET," S=",S
+#endif
+
+    if( dabs(DET) .lt. 1.0D-30)  then
+        RC = -1
+        AOBR(:,:)=0.0D00
+        return
+    endif
+
+    S=S/DET
+    AOBR(1,1)=A(2,2)*S
+    AOBR(1,2)=-A(1,2)*S
+    AOBR(2,1)=-A(2,1)*S
+    AOBR(2,2)=A(1,1)*S
+
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv2 : AOBR:"
+    write(ERROR_UNIT,'(2G20.12)') ((AOBR(i,j),j=1,2),i=1,2)
+#endif
+
+    return
+end
\ No newline at end of file

zve/math/nc_minv3.f95

+!>  @file nc_minv3.f95
+!>  @date Jun 15, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы 3x3
+!> @details  Подпрограмма выполняет обращение матрицы 3x3
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]    MA(3,3)    R*8 - обращаемая матрица A
+!> @param[out]   AOBR(3,3) R*8 - искомая обратная матрица
+!> @param[out]   RC        I*4 - код возврата
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+SUBROUTINE nc_minv3 ( MA, AOBR, RC )
+
+    use adaps_m
+
+    implicit none
+
+    real(8), intent(in)   :: MA(3,3)   !< обращаемая матрица
+    real(8), intent(out)  :: AOBR(3,3) !< обратная матрица
+    integer, intent(out)  :: RC        !< код возврата
+
+    real(8) :: nc_mdet3     ! функция поиска определителя матрицы 3x3
+
+    real(8) :: A(3,3)       ! Нормированная матрица
+    real(8) :: B(3,3)       ! Матрица алгебраических дополнений транспонированная
+    real(8) :: DET          ! определитель матрицы 3x3
+    real(8) :: S            ! множитель
+    
+    integer :: i,j          ! переменные цикла
+
+    RC = 0
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3 : MA:"
+    do i=1,3
+        write(ERROR_UNIT,'(3G20.12)') (MA(i,j),j=1,3)
+    enddo
+#endif
+
+    S=0.0D00
+    do i=1,3
+        do j=1,3
+            IF(DABS(MA(I,J)).gt.0.1D-50) S=S+DLOG(DABS(MA(i,j)))
+        enddo
+    enddo
+
+    S=S/9
+    S=1.0D00/DEXP(S)
+    
+    A(:,:)=MA(:,:)*S
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3 : A:"
+    do i=1,3
+        write(ERROR_UNIT,'(3G20.12)') (A(i,j),j=1,3)
+    enddo
+#endif
+
+
+    ! Находим алгебраические дополнения
+    B(1,1) =   A(2,2)*A(3,3)-A(2,3)*A(3,2)
+    B(2,1) = - A(2,1)*A(3,3)+A(2,3)*A(3,1)
+    B(3,1) =   A(2,1)*A(3,2)-A(2,2)*A(3,1)
+
+    B(1,2) = - A(1,2)*A(3,3)+A(1,3)*A(3,2)
+    B(2,2) =   A(1,1)*A(3,3)-A(1,3)*A(3,1)
+    B(3,2) = - A(1,1)*A(3,2)+A(1,2)*A(3,1)
+
+    B(1,3) =   A(1,2)*A(2,3)-A(1,3)*A(2,2)
+    B(2,3) = - A(1,1)*A(2,3)+A(1,3)*A(2,1)
+    B(3,3) =   A(1,1)*A(2,2)-A(1,2)*A(2,1)
+
+    ! Находим определитель 
+    DET = nc_mdet3(A)
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3 : DET=",DET," S=",S
+#endif
+
+    if( dabs(DET) .lt. 1.0D-30)  then
+        RC = -1
+        AOBR(:,:)=0.0D00
+        return
+    endif
+
+    S=S/DET
+    AOBR(:,:) = B(:,:)*S
+
+    return
+end
\ No newline at end of file

zve/math/nc_minv3f.f95

+!>  @file nc_minv3f.f95
+!>  @date Jun 25, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы 3x3 c использованием формулы Фробениуса
+!> @details  Подпрограмма выполняет обращение матрицы 3x3 c использованием
+!> формулы фробениуса
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]    M(3,3)    R*8 - обращаемая матрица M (туда же списывается результат обращения)
+!> @param[out]   RC        I*4 - код возврата
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+SUBROUTINE nc_minv3f ( MA, MOBR,RC )
+
+    use adaps_m
+
+    implicit none
+
+    real(8), intent(in)     :: MA(3,3)   !< обращаемая матрица
+    real(8), intent(out)    :: MOBR(3,3) !< обратная матрица
+    integer, intent(out)    :: RC        !< код возврата
+
+    real(8), parameter :: errtol = 1.0D-08 ! точность вычисления обратной матрицы
+
+    real(8) :: M(3,3)        ! нормированная обращаемая матрица
+    real(8) :: A1(1,1)       ! левый верхний блок матрицы M
+    real(8) :: B1(1,2)       ! правый верхний блок матрицы M
+    real(8) :: C1(2,1)       ! левый нижний блок матрицы M
+    real(8) :: D1(2,2)       ! правый нижний блок матрицы M
+    real(8) :: A2(2,2)       ! левый верхний блок матрицы M
+    real(8) :: B2(2,1)       ! правый верхний блок матрицы M
+    real(8) :: C2(1,2)       ! левый нижний блок матрицы M
+    real(8) :: D2(1,1)       ! правый нижний блок матрицы M
+    real(8) :: AOBR1(1,1)    ! обратная матрица для левого верхнего блока
+    real(8) :: AOBR2(2,2)    ! обратная матрица для левого верхнего блока
+    real(8) :: H1(2,2)       ! вспомогательная матрица
+    real(8) :: H2(1,1)       ! вспомогательная матрица
+    real(8) :: HOBR1(2,2)    ! обратная матрица к вспомогательной матрицы
+    real(8) :: HOBR2(1,1)    ! обратная матрица к вспомогательной матрицы
+    real(8) :: W1(1,1)       ! рабочая матрица
+    real(8) :: W2(2,2)       ! рабочая матрица
+    real(8) :: W12(1,2)      ! рабочая матрица
+    real(8) :: W21(2,1)      ! рабочая матрица
+    real(8) :: E(3,3)        ! единичная матрица
+    real(8) :: MEM(3,3)      ! обратная матрица запомненная
+
+    real(8) :: maxerr1       ! максимальное значение ошибки
+    real(8) :: maxerr2       ! максимальное значение ошибки на второй итерации
+    real(8) :: S             ! множитель
+    
+    real(8) :: R2(2,2)       ! рабочая матрица
+    real(8) :: R1(1,1)       ! рабочая матрица
+    
+    integer :: i,j          ! переменные цикла
+#ifdef DEBUG        
+    real(8) :: E2(2,2)       ! единичная матрица
+#endif
+    
+    RC = 0
+    MOBR(:,:) = 0.0D00
+
+    S=0.0D00
+    do i=1,3
+        do j=1,3
+            IF(DABS(MA(I,J)).gt.0.1D-50) S=S+DLOG(DABS(MA(i,j)))
+        enddo
+    enddo
+
+    S=S/9
+    S=1.0D00/DEXP(S)
+
+    M(:,:)=MA(:,:)*S
+
+    ! выполняем присвоения
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3f : M:"
+    write(ERROR_UNIT,'(3G20.12)') ((M(i,j),j=1,3),i=1,3)
+#endif
+    A1(1,1)=M(1,1)
+    do i=1,2
+        B1(1,i)=M(1,i+1)
+        C1(i,1)=M(i+1,1)
+    enddo
+    do i=1,2
+        do j=1,2
+            D1(i,j)=M(i+1,j+1)
+        enddo
+    enddo
+    
+    ! 1 вариант
+    !           A1(1x1) B1(1x2)
+    ! M(3x3) = 
+    !           C1(2x1) D1(2x2)
+
+    AOBR1(1,1)=1.0D00/A1(1,1)
+
+    ! W2(2x2) = C1(2x1)*A(1x1)^(-1)*B1(1x2)
+    call nc_mprm3(2,1,1,2,C1,AOBR1,B1,W2)
+    
+    ! H2(2x2) = D1(2x2)-W2(2x2)
+    do i=1,2
+        do j=1,2
+            H1(i,j)=D1(i,j)-W2(i,j)
+        enddo
+    enddo
+    
+    ! HOBR1(2x2)=H1(2x2)^(-1)
+    call nc_minv2(H1,HOBR1,RC)
+
+#ifdef DEBUG
+    ! Проверка
+    call nc_mprm(2,2,2,H1,HOBR1,E2)
+
+    write(ERROR_UNIT,*) "nc_minv3f : HOBR1:"
+    write(ERROR_UNIT,'(2G20.12)') ((HOBR1(i,j),j=1,2),i=1,2)
+
+    write(ERROR_UNIT,*) "nc_minv3f : E2:"
+    write(ERROR_UNIT,'(2G20.12)') ((E2(i,j),j=1,2),i=1,2)
+#endif
+    if(RC.ne.0) return
+    
+    ! W1(1x1)=B(1x2)*HOBR1(2x2)*C(2x1)
+    call nc_mprm3(1,2,2,1,B1,HOBR1,C1,W1)
+    
+    ! R1(1x1)=AOBR1(1x1)*W1(1,1)*AOBR1(1,1)
+    R1(1,1)=AOBR1(1,1)*W1(1,1)*AOBR1(1,1)
+    
+    ! W12(1x2)=AOBR1(1x1)*B1(1x2)*HOBR1(2,2)
+    call nc_mprm3(1,1,2,2,AOBR1,B1,HOBR1,W12)
+    
+    ! W21(2x1)=HOBR1(2x2)*C(2x1)*AOBR(1,1)
+    call nc_mprm3(2,2,1,1,HOBR1,C1,AOBR1,W21)
+    
+    MOBR(1,1)=AOBR1(1,1)+R1(1,1)
+    do i=1,2
+        MOBR(1,i+1)=-W12(1,i)
+        MOBR(i+1,1)=-W21(i,1)
+    enddo    
+    do i=1,2
+        do j=1,2
+            MOBR(i+1,j+1)=HOBR1(i,j)
+        enddo
+    enddo
+    
+    ! Проверка:
+    call nc_mprm(3,3,3,MOBR,M,E)
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3f : E:"
+    write(ERROR_UNIT,'(3G20.12)') ((E(i,j),j=1,3),i=1,3)
+#endif
+
+    ! находим максимальную ошибку
+    maxerr1 = 0.0D00
+    do i=1,3
+        do j=1,3
+            if(i.ne.j) then
+                if(dabs(E(i,j)).gt.maxerr1) maxerr1 = dabs(E(i,j))
+            else
+                if(dabs(E(i,j) - 1.0D00).gt.maxerr1) maxerr1 = dabs(E(i,j)-1.0D00)
+            endif
+        enddo
+    enddo
+    if(maxerr1.lt.errtol) then
+        MOBR(:,:) = MOBR(:,:)*S
+        return
+    endif
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3f : maxerr1:",maxerr1," --> new attempt"
+#endif
+
+    MEM(:,:)=MOBR(:,:)
+
+    ! 2 вариант
+    !           A2(2x2) B2(2x1)
+    ! M(3x3) = 
+    !           C2(1x2) D2(1x1)
+
+    ! AOBR2(2x2) = A2(2x2)^(-1)
+
+
+    do i=1,2
+        do j=1,2 
+            A2(i,j)=M(i,j)
+        enddo
+    enddo
+    do i=1,2        
+        B2(i,1)=M(i,3)
+        C2(1,i)=M(3,i)
+    enddo    
+    D2(1,1)=M(3,3)
+
+    call nc_minv2(A2,AOBR2,RC)
+
+#ifdef DEBUG
+    ! Проверка
+    call nc_mprm(2,2,2,A2,AOBR2,E2)
+
+    write(ERROR_UNIT,*) "nc_minv3f : AOBR1:"
+    write(ERROR_UNIT,'(2G20.12)') ((AOBR2(i,j),j=1,2),i=1,2)
+
+    write(ERROR_UNIT,*) "nc_minv3f : E2:"
+    write(ERROR_UNIT,'(2G20.12)') ((E2(i,j),j=1,2),i=1,2)
+#endif
+    if(RC.ne.0) return
+
+    ! W1(1x1) = C2(1x2)*A(1x1)^(-1)*B2(2x1)
+    call nc_mprm3(1,2,2,1,C2,AOBR2,B2,W1)
+    
+    ! H2(1x1) = D1(1x1)-W1(1x1)
+    H2(1,1)=D2(1,1)-W1(1,1)
+
+    ! HOBR2(1x1)=1/H2(1x1)
+    HOBR2(1,1)=1.0D00/H2(1,1)
+
+    ! W2(2x2)=B2(2x1)*HOBR2(1x1)*C2(1x2)
+    call nc_mprm3(2,1,1,2,B2,HOBR2,C2,W2)
+
+    ! R2(2x2)=AOBR2(2x2)*W2(2,2)*AOBR2(2,2)
+    call nc_mprm3(2,2,2,2,AOBR2,W2,AOBR2,R2)
+
+    ! W21(2x1)=AOBR2(2x2)*B2(2x1)*HOBR2(1,1)
+    call nc_mprm3(2,2,1,1,AOBR2,B2,HOBR2,W21)
+
+    ! W12(1x2)=HOBR2(1x1)*C2(1x2)*AOBR2(2,2)
+    call nc_mprm3(1,1,2,2,HOBR2,C2,AOBR2,W12)
+
+    do i=1,2
+        do j=1,2
+            MOBR(i,j)=AOBR2(i,j)+R2(i,j)
+        enddo
+    enddo        
+    do i=1,2
+        MOBR(i,3)=-W21(i,1)
+        MOBR(3,i)=-W12(1,i)
+    enddo
+    MOBR(3,3)=HOBR2(1,1)    
+
+    ! Проверка:
+    call nc_mprm(3,3,3,MOBR,M,E)
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3f : 2 attempt E:"
+    write(ERROR_UNIT,'(3G20.12)') ((E(i,j),j=1,3),i=1,3)
+#endif
+    ! находим максимальную ошибку
+    maxerr2 = 0.0D00
+    do i=1,3
+        do j=1,3
+            if(i.ne.j) then
+                if(dabs(E(i,j)).gt.maxerr2) maxerr2 = dabs(E(i,j))
+            else
+                if(dabs(E(i,j) - 1.0D00).gt.maxerr2) maxerr2 = dabs(E(i,j)-1.0D00)
+            endif
+        enddo
+    enddo
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv3f : maxerr2:",maxerr2
+#endif
+
+    if(maxerr1.lt.maxerr2) then
+        MOBR(:,:)=MEM(:,:)
+    endif    
+
+    MOBR(:,:) = MOBR(:,:)*S
+    
+    return
+end
\ No newline at end of file

zve/math/nc_minv6.f95

+!>  @file nc_minv6.f95
+!>  @date Jun 25, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы 6x6
+!> @details  Подпрограмма выполняет обращение матрицы 6x6 c использованием
+!> формулы фробениуса
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]    M(6,6)    R*8 - обращаемая матрица M (туда же списывается результат обращения)
+!> @param[out]   RC        I*4 - код возврата
+!>
+!> @retval  0   нориальное завершение
+!> @retval -1   определитель матрицы равен 0
+!>
+SUBROUTINE nc_minv6 ( M, RC )
+
+    use adaps_m
+
+    implicit none
+    
+    real(8), parameter :: errtol = 1.0D-08 ! Заданная точность 
+
+    real(8), intent(inout)  :: M(6,6)    !< обращаемая матрица
+    integer, intent(out)    :: RC        !< код возврата
+
+    real(8) :: A(3,3)       ! левый верхний блок матрицы M
+    real(8) :: B(3,3)       ! правый верхний блок матрицы M
+    real(8) :: C(3,3)       ! левый нижний блок матрицы M
+    real(8) :: D(3,3)       ! правый нижний блок матрицы M
+    real(8) :: AOBR(3,3)    ! обратная матрица для левого верхнего блока
+    real(8) :: H(3,3)       ! вспомогательная матрица
+    real(8) :: HOBR(3,3)    ! обратная матрица к вспомогательной матрицы
+    real(8) :: WOBR(3,3)    ! обратная матрица к вспомогательной матрицы
+    real(8) :: W(3,3)       ! рабочая матрица
+    real(8) :: R(3,3)       ! рабочая матрица
+    real(8) :: T(3,3)       ! рабочая матрица
+    real(8) :: E(6,6)       ! единичная матрица
+    real(8) :: O(6,6)       ! оригинальная матрица
+    real(8) :: maxerr1      ! максимальное значение ошибки
+    real(8) :: maxerr2      ! максимальное значение ошибки на второй итерации
+    
+    integer :: i,j          ! переменные цикла
+    
+    RC = 0
+
+    ! выполняем присвоения
+    do i=1,3
+        do j=1,3
+            A(i,j)=M(i,j)
+            B(i,j)=M(i,j+3)
+            C(i,j)=M(i+3,j)
+            D(i,j)=M(i+3,j+3)
+        enddo
+    enddo
+    
+    O(:,:) = M(:,:)
+    M(:,:) = 0.0D00
+    
+    call nc_minv3(A,AOBR,RC)
+        
+    call nc_mprm(3,3,3,A,AOBR,W)
+#ifdef DEBUG
+        write(ERROR_UNIT,*) "nc_minv6 : nc_minv3 : A*AOBR = ONE MATRIX :"
+        write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+    
+    maxerr1 = 0.0D00
+    do i=1,3
+        do j=1,3
+            if(i.ne.j) then
+                if(dabs(W(i,j)).gt.maxerr1) maxerr1 = dabs(W(i,j))
+            else
+                if(dabs(W(i,j) - 1.0D00).gt.maxerr1) maxerr1 = dabs(W(i,j)-1.0D00)
+            endif
+        enddo
+    enddo
+
+    if(maxerr1.gt.errtol) then
+        WOBR(:,:)=AOBR(:,:)
+        call nc_minv3f(A,AOBR,RC)
+        call nc_mprm(3,3,3,A,AOBR,W)
+#ifdef DEBUG
+        write(ERROR_UNIT,*) "nc_minv6 : nc_minv3f : A*AOBR = ONE MATRIX :"
+        write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+        if(RC.ne.0) then
+            write(ERROR_UNIT,*) "nc_minv6 : DEBUG : inversion A : RC=",RC
+            return
+        endif    
+        maxerr2 = 0.0D00
+        do i=1,3
+            do j=1,3
+                if(i.ne.j) then
+                    if(dabs(W(i,j)).gt.maxerr2) maxerr2 = dabs(W(i,j))
+                else
+                    if(dabs(W(i,j) - 1.0D00).gt.maxerr2) maxerr2 = dabs(W(i,j)-1.0D00)
+                endif
+            enddo
+        enddo
+        if(maxerr2.gt.maxerr1) then
+#ifdef DEBUG
+            write(ERROR_UNIT,*) "nc_minv6 : nc_minv3 solution AOBR:"
+            write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+            AOBR(:,:)=WOBR(:,:)
+        endif
+    endif    
+
+#ifdef DEBUG
+    ! Проверка:
+    call nc_mprm(3,3,3,A,AOBR,W)
+    write(ERROR_UNIT,*) "nc_minv6 : A*AOBR = ONE MATRIX :"
+    write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+
+    call nc_mprm3(3,3,3,3,C,AOBR,B,W)
+    do i=1,3
+        do j=1,3
+            H(i,j)=D(i,j)-W(i,j)
+        enddo
+    enddo
+
+    call nc_minv3(H,HOBR,RC)
+
+    call nc_mprm(3,3,3,H,HOBR,W)
+#ifdef DEBUG
+        write(ERROR_UNIT,*) "nc_minv6 : nc_minv3 : H*HOBR = ONE MATRIX :"
+        write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+    maxerr1 = 0.0D00
+    do i=1,3
+        do j=1,3
+            if(i.ne.j) then
+                if(dabs(W(i,j)).gt.maxerr1) maxerr1 = dabs(W(i,j))
+            else
+                if(dabs(W(i,j) - 1.0D00).gt.maxerr1) maxerr1 = dabs(W(i,j)-1.0D00)
+            endif
+        enddo
+    enddo
+    if(maxerr1.gt.errtol) then
+        WOBR(:,:)=HOBR(:,:)
+        call nc_minv3f(H,HOBR,RC)
+        call nc_mprm(3,3,3,H,HOBR,W)
+#ifdef DEBUG
+        write(ERROR_UNIT,*) "nc_minv6 : nc_minv3f : H*HOBR = ONE MATRIX :"
+        write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+        if(RC.ne.0) then
+            write(ERROR_UNIT,*) "nc_minv6 : DEBUG : inversion H : RC=",RC
+            return
+        endif    
+        maxerr2 = 0.0D00
+        do i=1,3
+            do j=1,3
+                if(i.ne.j) then
+                    if(dabs(W(i,j)).gt.maxerr2) maxerr2 = dabs(W(i,j))
+                else
+                    if(dabs(W(i,j) - 1.0D00).gt.maxerr2) maxerr2 = dabs(W(i,j)-1.0D00)
+                endif
+            enddo
+        enddo
+        if(maxerr2.gt.maxerr1) then
+#ifdef DEBUG
+            write(ERROR_UNIT,*) "nc_minv6 : nc_minv3 solution HOBR:"
+            write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+            HOBR(:,:)=WOBR(:,:)
+        endif
+    endif    
+
+#ifdef DEBUG
+    ! Проверка
+    call nc_mprm(3,3,3,H,HOBR,W)
+    write(ERROR_UNIT,*) "nc_minv6 : H*HOBR = ONE MATRIX :"
+    write(ERROR_UNIT,'(3G20.12)') ((W(i,j),j=1,3),i=1,3)
+#endif
+    
+    call nc_mprm3(3,3,3,3,B,HOBR,C,W)
+    call nc_mprm3(3,3,3,3,AOBR,W,AOBR,R)
+    call nc_mprm3(3,3,3,3,AOBR,B,HOBR,W)
+    call nc_mprm3(3,3,3,3,HOBR,C,AOBR,T)
+    
+    do i=1,3
+        do j=1,3
+            M(i,j)=AOBR(i,j)+R(i,j)
+            M(i,j+3)=-W(i,j)
+            M(i+3,j)=-T(i,j)
+            M(i+3,j+3)=HOBR(i,j)
+        enddo
+    enddo
+
+    ! Проверка    
+    call nc_mprm(6,6,6,O,M,E)
+#ifdef DEBUG
+    write(ERROR_UNIT,*) "nc_minv6 : rezult E:"
+    write(ERROR_UNIT,'(6G20.12)') ((E(i,j),j=1,6),i=1,6)
+#endif
+    
+    return
+end
\ No newline at end of file

zve/math/nc_minvn.f95

+!>  @file nc_minvn.f95
+!>  @date Jun 18, 2020
+!>
+!> @brief  Подпрограмма обращения матрицы c нормированием
+!> @details  Подпрограмма выполняет обращение матрицы с нормированием
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]      N       I*4 - размерность матрицы А
+!> @param[inout]   A(N,N)  R*8 - обращаемая матрица A (туда же списывается результат обращения)
+!> @param[out]     DET     R*4 - определитель матрицы
+!> @param[inout]   LK(N)   I*4 - вспомогательный целочисленный массив
+!> @param[inout]   KL(N)   I*4 - вспомогательный целочисленный массив
+!>
+SUBROUTINE nc_minvn(N,DET,A,LK,KL)
+
+    use adaps_m
+
+    implicit none
+
+    real(8), intent(inout) :: A(N,N)  !< обращаемая матрица
+    integer, intent(in)    :: N       !< размерность матрицы
+    real(8), intent(out)   :: DET     !< определитель матрицы
+    integer, intent(inout) :: LK(N)   !< вспомогательный целочисленный массив
+    integer, intent(inout) :: KL(N)   !< вспомогательный целочисленный массив
+    
+    real(8) :: S                      ! нормирующий множитель
+    integer :: i,j                    ! переменные цикла
+
+    if(N.le.0) then
+        write(ERROR_UNIT,*) "nc_minvn : DEBUG: control dimension"
+    endif
+
+    DET = 0.0D00
+    S=0.0D00
+    do i=1,N
+        do j=1,N
+            IF(DABS(A(I,J)).gt.0.1D-50) S=S+DLOG(DABS(A(i,j)))
+        enddo
+    enddo
+
+    S=S/(N*N)
+    S=1.0D00/DEXP(S)
+    
+    do i=1,N
+        do j=1,N
+            A(i,j)=A(i,j)*S
+        enddo
+    enddo
+
+    CALL nc_dminv(A,N,DET,LK,KL)
+    
+    do i=1,N
+        do j=1,N
+            A(i,j)=A(i,j)*S
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_mmov.f95

+!>  @file nc_mmov.f95
+!>  @date Jun 21, 2020
+!>
+!> @brief  Копироавние содержимого матрицы 
+!> @details  Подпрограмма выполняет копирование содержимого матрицы А в B
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]     N   I*4 - количество строк в матрице
+!> @param[in]     M   I*4 - количество столбцов в матрице
+!> @param[in]     A   R*8 - матрица A(NxM)
+!> @param[out]    B   R*8 - матрица A(NxM)
+!>
+SUBROUTINE nc_mmov(N,M,A,B)
+
+    use adaps_m
+
+    implicit none
+
+    integer, intent(in)      :: N       !< размерность матрицы
+    integer, intent(in)      :: M       !< размерность матрицы
+    real(8), intent(in)      :: A(N,M)  !< копируемая матрица A
+    real(8), intent(out)     :: B(N,M)  !< матрица B
+
+    integer :: I,J                      ! переменные цикла
+
+    do I=1,N
+        do J=1,M
+            B(I,J) = A(I,J)
+        enddo
+    enddo        
+    RETURN
+END
\ No newline at end of file

zve/math/nc_mone.f95

+!>  @file nc_mone.f95
+!>  @date Jan 16, 2022
+!>
+!> @brief  Единичная матрица
+!> @details  Подпрограмма формирует единичную матрицу
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]     N   I*4 - количество строк в матрице
+!> @param[inout]  A   R*8 - матрица A(NxM)
+!>
+SUBROUTINE nc_mone(N,A)
+
+    use adaps_m
+
+    implicit none
+
+    integer, intent(in)      :: N       !< размерность матрицы
+    real(8), intent(inout)   :: A(N,N)  !< единичная матрица
+
+    integer :: I,J                      ! переменные цикла
+
+    do i=1,N
+        do j=1,N
+            if(i.ne.j) then
+                A(i,j) = 0.0D00
+            else
+                A(i,j) = 1.0D00
+            endif
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_mprm.f95

+!>  @file nc_mprm.f95
+!>  @date Jun 18, 2020
+!>
+!> @brief  Подпрограмма перемножения матриц с нормированием
+!> @details  Подпрограмма выполняет перемножение двух матриц с нормированием
+!> R(NxL)=A(NxM)*B(MxL)
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   N       I*4 - количество строк в матрице А
+!> @param[in]   M       I*4 - количество строк в матрице B и столбцов в матрице A
+!> @param[in]   L       I*4 - количество столбцов в матрице B
+!> @param[in]  A(N,M)   R*8 - матрица A
+!> @param[in]  B(M,L)   R*8 - матрица B
+!> @param[out] R(N,L)   R*8 - произведение матриц A и B
+!>
+
+SUBROUTINE nc_mprm(N,M,L,A,B,R)
+
+    use adaps_m
+    
+    IMPLICIT NONE
+    
+    integer, intent(in)  :: N       !< количество строк в матрице А
+    integer, intent(in)  :: M       !< количество строк в матрице B и столбцов в матрице A
+    integer, intent(in)  :: L       !< количество столбцов в матрице B
+    REAL(8), intent(in)  :: A(N,M)  !< матрица A
+    REAL(8), intent(in)  :: B(M,L)  !< матрица B
+    REAL(8), intent(out) :: R(N,L)  !< матрица R = A*B
+
+    REAL(8) :: SN1     ! среднее значение логарифмов элементов матрицы А
+    REAL(8) :: SN2     ! среднее значение логарифмов элементов матрицы B
+    REAL(8) :: AM1     ! нормирующий множитель матрицы A
+    REAL(8) :: AM2     ! нормирующий множитель матрицы B
+    REAL(8) :: AM      ! нормирующий множитель матрицы R
+    REAL(8) :: PR      ! рабочая переменная 
+    
+    integer :: i,j,k   ! переменные цикла
+
+    if(N.le.0.or.M.le.0.or.L.le.0) then
+        write(ERROR_UNIT,*) "nc_nvmprm : DEBUG: control dimension"
+    endif
+
+    SN1=0.0D00
+    do i=1,N
+        do j=1,M
+            IF(DABS(A(i,j)).gt.0.1D-50) SN1=SN1+DLOG(DABS(A(i,j)))
+        enddo
+    enddo
+    
+    SN2=0.0D00
+    do i=1,M
+        do j=1,L
+            IF(DABS(B(i,j)).gt.0.1D-50) SN2=SN2+DLOG(DABS(B(i,j)))
+        enddo
+    enddo        
+
+    SN1=SN1/(N*M)
+    SN2=SN2/(M*L)
+    AM1=1.0D00/DEXP(SN1)
+    AM2=1.0D00/DEXP(SN2)
+    AM=DEXP(SN1+SN2)
+    
+    do i=1,N
+        do k=1,L
+            PR=0.0D00
+            do j=1,M
+                PR=PR+(AM1*A(i,j))*(AM2*B(j,k))
+            enddo
+            R(i,k)=PR*AM
+        enddo
+    enddo
+    RETURN
+END

zve/math/nc_mprm3.f95

+!>  @file nc_mprm3.f95
+!>  @date Jun 18, 2020
+!>
+!> @brief  Подпрограмма перемножения трех матриц с нормированием
+!> @details  Подпрограмма выполняет перемножение трех матриц с нормированием
+!> R(NxK)=A(NxM)*B(MxL)*С(L,K)
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   N       I*4 - количество строк в матрице А
+!> @param[in]   M       I*4 - количество строк в матрице B и столбцов в матрице A
+!> @param[in]   L       I*4 - количество строк в матрице С и столбцов в матрице B
+!> @param[in]   K       I*4 - количество столбцов в матрице C
+!> @param[in]  A(N,M)   R*8 - матрица A
+!> @param[in]  B(M,L)   R*8 - матрица B
+!> @param[in]  C(L,K)   R*8 - матрица C
+!> @param[out] R(N,K)   R*8 - произведение матриц A и B и С
+!>
+
+SUBROUTINE nc_mprm3(N,M,L,K,A,B,C,R)
+
+    use adaps_m
+    
+    IMPLICIT NONE
+    
+    integer, intent(in)  :: N       !< количество строк в матрице А
+    integer, intent(in)  :: M       !< количество строк в матрице B и столбцов в матрице A
+    integer, intent(in)  :: L       !< количество строк в матрице C и столбцов в матрице B
+    integer, intent(in)  :: K       !< количество столбцов в матрице C
+    REAL(8), intent(in)  :: A(N,M)  !< матрица A
+    REAL(8), intent(in)  :: B(M,L)  !< матрица B
+    REAL(8), intent(in)  :: C(L,K)  !< матрица C
+    REAL(8), intent(out) :: R(N,K)  !< матрица R = A*B*C
+
+    REAL(8) :: SN1       ! среднее значение логарифмов элементов матрицы А
+    REAL(8) :: SN2       ! среднее значение логарифмов элементов матрицы B
+    REAL(8) :: SN3       ! среднее значение логарифмов элементов матрицы C
+    REAL(8) :: AM1       ! нормирующий множитель матрицы A
+    REAL(8) :: AM2       ! нормирующий множитель матрицы B
+    REAL(8) :: AM3       ! нормирующий множитель матрицы C
+    REAL(8) :: AM        ! нормирующий множитель матрицы R
+    REAL(8) :: PR        ! рабочая переменная 
+    
+    integer :: i,j       ! переменные цикла
+    integer :: i1,i2,i3  ! переменные цикла
+
+    if(N.le.0.or.M.le.0.or.L.le.0.or.K.le.0) then
+        write(ERROR_UNIT,*) "nc_nvmprm3 : DEBUG: control dimension"
+    endif
+
+    SN1=0.0D00
+    do i=1,N
+        do j=1,M
+            if(DABS(A(i,j)).gt.0.1D-50) SN1=SN1+DLOG(DABS(A(i,j)))
+        enddo
+    enddo
+
+    SN2=0.0D00
+    do i=1,M
+        do j=1,L
+            if(DABS(B(i,j)).gt.0.1D-50) SN2=SN2+DLOG(DABS(B(i,j)))
+        enddo
+    enddo
+
+    SN3=0.0D00
+    do i=1,L
+        do j=1,K
+            if(DABS(C(i,j)).gt.0.1D-50) SN3=SN3+DLOG(DABS(C(i,j)))
+        enddo
+    enddo
+
+    SN1=SN1/(N*M)
+    SN2=SN2/(M*L)
+    SN3=SN3/(L*K)
+    AM1=1.0D00/DEXP(SN1)
+    AM2=1.0D00/DEXP(SN2)
+    AM3=1.0D00/DEXP(SN3)
+    AM=DEXP(SN1+SN2+SN3)
+
+    CALL nc_mnull(N,K,R)
+
+    do i1=1,N
+        do i2=1,L
+            PR=0.0D00
+            do i3=1,M
+                PR=PR+(AM1*A(i1,i3))*(AM2*B(i3,i2))
+            enddo    
+            do i3=1,K
+                R(i1,i3)=R(i1,i3)+PR*AM3*C(i2,i3)
+            enddo
+        enddo
+    enddo
+
+    do i1=1,N
+        do i2=1,K
+            R(i1,i2)=R(i1,i2)*AM
+        enddo
+    enddo
+
+    RETURN
+END

zve/math/nc_msim.f95

+!>  @file nc_msim.f95
+!>  @date Jun 21, 2020
+!>
+!> @brief  Решение системы линейных алгебраических уравнений
+!> @details  Подпрограмма решает систему линейных алгебраических уравнений АX=B
+!> методом исключения.
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[inout]   A(NxN)   R*8 - матрица коэффициентов порядка N x N, расположенная по столбцам.
+!>                                В ходе вычисления изменяется.
+!> @param[inout]   B(N)     R*8 - вектор свободных членов системы длины N,
+!>                                На выходе заменяется вектором решения X.
+!> @param[in]      N        I*4 - количество уравнений, N > 1.
+!> @param[out]     KS       I*4 - признак решения принимает значения:
+!>                                0 - для нормального решения,
+!>                                1 - для особенной системы уравнений.
+!>
+!> 
+!> @note
+!>  Maтрица А должна быть общего вида. Если матрица вырожденная, то решение не имеет смысла.
+!>  Решение может быть получено также путем использования программы обращения матрицы (MINV)
+!>  с последующим расчётом произведения матриц.
+!>
+!> Описание метода:
+!>
+!>  Peшение производится методом исключения с выбором главного элемента. При необходимости каждый шаг 
+!>  исключения содержит перестановку строк, что бы избежать на нуль или очень маленький элемент.
+!>  k-я неизвестная получается при прямом ходе решения за N шагов. Остальные неизвестные получаются
+!>  с помощью последовательных подстанов ок при обратном ходе решения.
+!>
+!>  Решение размещается в векторе В: первая неизвестная в В(1), вторая - в В(2), N-я неизвестная - в В(N).
+!>
+!>  Если главный элемент не превышает допуска TOL=0.0, то матрица принимается за вырожденную и MS предполагается равным 1.
+!>
+!>  Допуск TOL можно изменить, явно изменяя его значение.
+!>
+subroutine nc_msim(A,B,N,KS)
+    use adaps_m
+
+    implicit none
+
+    integer, intent(in)      :: N       !< количество уравнений
+    integer, intent(out)     :: KS      !< признак решения
+    real(8), intent(inout)   :: A(N)    !< матрица коэффициентов порядка N x N, расположенная по столбцам.
+    real(8), intent(inout)   :: B(N)    !< вектор свободных членов системы длины N,
+
+    REAL(8) ::  TOL            ! Допустимый порог
+    REAL(8) ::  BIGA           ! Наибольшее значение матрицы
+    REAL(8) ::  SAV            ! сохраненное значение
+    
+    integer :: i,j,k           ! переменные цикла
+    integer :: ix,jx           ! переменные цикла
+    integer :: imax,ny,ij      ! рабочая переменная
+    integer :: jj,ixjx,jjx     ! рабочая переменная
+    integer :: jy,ixj,iqs      ! рабочая переменная
+    integer :: it,i1,i2        ! рабочая переменная
+    integer :: ia,ib,ic        ! рабочая переменная
+
+    TOL=0.0D00
+    KS=0
+    imax=0
+    jj=-N
+    do j=1,N
+        jy=j+1
+        jj=jj+N+1
+        BIGA=0.0D00
+        it=jj-j
+        do i=j,N
+            ij=it+i
+            IF( (DABS(BIGA)-DABS(A(ij)) ).lt.0.0D00) then
+                BIGA=A(ij)
+                imax=i
+            endif
+        enddo
+        if((DABS(BIGA)-TOL).le.0.0D00) then
+            KS=1
+            return
+        endif
+        i1=j+N*(j-2)
+        it=imax-j
+        do k=j,N
+            i1=i1+N
+            i2=i1+it
+            SAV=A(i1)
+            A(i1)=A(i2)
+            A(i2)=SAV
+            A(i1)=A(i1)/BIGA
+        enddo
+        SAV=B(imax)
+        B(imax)=B(j)
+        B(j)=SAV/BIGA
+        if((j-N).eq.0) goto 70
+        iqs=N*(j-1)
+        do ix=jy,N
+            ixj=iqs+ix
+            it=j-ix
+            do jx=jy,N
+                ixjx=N*(jx-1)+IX
+                jjx=IXJX+it
+                A(ixjx)=A(ixjx)-(A(ixj)*A(jjx))
+            enddo
+            B(ix)=B(ix)-(B(j)*A(ixj))
+        enddo
+    enddo
+
+70  ny=N-1
+    it=N*N
+    do j=1,NY
+        ia=it-j
+        ib=N-j
+        ic=N
+        do k=1,j
+            B(ib)=B(ib)-A(ia)*B(ic)
+            ia=ia-N
+            ic=ic-1
+        enddo
+    enddo
+    RETURN
+END
+

zve/math/nc_ndtr.f95

+!>  @file nc_ndtr.f95
+!>  @date Jun 18, 2020
+!>
+!> @brief  Вычисляется значение функции Лапласа
+!> @details  Подпрограмма вычисления функции Лапласа
+!>
+!> Фopмaльныe пapaмeтpы: 
+!>
+!> @param[in]   X   R*8 - скалярная величина для которой вычисляется функция Лапласа
+!> @param[out]  P   R*8 - значение интеграла Лапласа
+!> @param[out]  P   R*8 - плотность вероятности ( значение подынтегральной функции)
+!>
+SUBROUTINE nc_ndtr(X,P,D)
+
+    implicit none
+
+    REAL(8), intent(in)  :: X  !< скалярная величина для которой вычисляется функция Лапласа
+    REAL(8), intent(out) :: P  !< значение интеграла Лапласа
+    REAL(8), intent(out) :: D  !< плотность вероятности ( значение подынтегральной функции)
+
+    REAL(8) :: AX   ! модуль скалярной величины Х
+    REAL(8) :: T    ! рабочая переменная 
+
+    AX=DABS(X)
+    T = 1.0D00/(1.0D00+0.2316419D00*AX)
+    D = 0.3989423D00*DEXP(-X*X/2.0D00)
+    P = 1.0D00 - D*T*((((1.330274D00*T - 1.821256D00)*T + 1.781478D00)*T - 0.3565638D00)*T + 0.3193815D00)
+    IF(X.lt.0.0D00) P = 1.0D00-P
+    RETURN
+END

zve/math/nvmPLANE.f

+! ------------------------------------------------------------------------------
+!
+!                           SUBROUTINE nvmPLANE
+!
+! ------------------------------------------------------------------------------
+!
+!  Подпрограмма определения плоскости по известным трем точкам
+!   a b c d переменные описывающие плоскость. 
+!   Уравнение плоскости :  ax + by + cz + d =ь0 
+!   искомые коэффициенты определяются при помощи нахождения определителя матрицы:
+!                                                    x  y  z   1
+!                                                    x1 y1 z1  1   =0
+!                                                    x2 y2 z2  1
+!                                                    x3 y3 z3  1
+!  Формальные параметры:
+!    X1(3) -    x1,y1,z1    - точка 1
+!    X2(3) -    x2,y2,z2    - точка 2
+!    X3(3) -    x3,y3,z3    - точка 3
+!    a,b,c,d     - коэффициенты в уравнении плоскости
+!
+! ------------------------------------------------------------------------------  
+
+      SUBROUTINE nvmPLANE  ( X1,X2,X3, a,b,c,d )
+      IMPLICIT NONE
+      real*8 X1(3),X2(3),X3(3)
+      REAL*8 a,b,c,d
+
+      REAL*8 z23,y23,x23,yz23,yz32,xz23,xz32,xy23,xy32
+
+      z23 = X2(3)-X3(3)  ! z23= z2-z3 
+      y23 = X2(2)-X3(2)  ! y23= y2-y3 
+      x23 = X2(1)-X3(1)  ! x23= x2-x3 
+
+      yz23 = X2(2)*X3(3) ! yz23= y2*z3 
+      yz32 = X3(2)*X2(3) ! yz32= y3*z2 
+      xz23 = X2(1)*X3(3) ! xz23= x2*z3
+      xz32 = X3(1)*X2(3) ! xz32= x3*z2 
+      xy23 = X2(1)*X3(2) ! xy23= x2*y3 
+      xy32 = X3(1)*X2(2) ! xy32= x3*y2 
+
+      a= X1(2)*z23 - X1(3)*y23 + yz23 - yz32 
+      b= -(X1(1)*z23 - X1(3)*x23 + xz23 - xz32)
+      c= X1(1)*y23 - X1(2)*x23 + xy23 - xy32 
+      d= -(X1(1)*(yz23 - yz32) - X1(2)*(xz23-xz32) + X1(3)*(xy23 -xy32))
+      RETURN
+      END

zve/math/nvmkpl.f

+      REAL*8 FUNCTION nvmKPL(EX,M)
+*---------------------------------------------------------------+
+* Вычисление эксцентрической аномалии по известным
+* средней аномалии и эксцентриситету
+*---------------------------------------------------------------+
+* Bxoдныe пapaмeтpы: !
+* EX - эксцентриситет
+* M - средняя аномалия
+* Функция возвращает значение эксцентрической аномалии (рад) [0,2PI]
+* Решение осуществляется методом последовательных приближений
+*---------------------------------------------------------------+
+      IMPLICIT NONE
+      REAL*8  EX,M,E0,EK
+      REAL*8  PREC/1.0D-12/
+
+*  Задаем точность и начальное приближение
+      E0=M
+*  Решение методом последовательных приближений
+1     EK=EX*DSIN(E0)+M
+      IF(DABS(EK-E0).LT.PREC) GOTO 99999
+      E0=EK
+      GOTO 1
+
+99999 nvmKPL=EK
+      return
+      end

zve/mod/.gitignore

+!Doxyfile
\ No newline at end of file

zve/mod/harmo.f

+       SUBROUTINE HARMO(NN,XB,RMU,RR,C,U,GRAD)
+       IMPLICIT REAL*8(A-H,O-Z)
+       DIMENSION XB(3),X(3),GRAD(3),C(1)
+       DATA Q0,Q1/0.D0,1.D0/,EPS/2.D-12/
+       RO=XB(1)*XB(1)+XB(2)*XB(2)
+       RV=DSQRT(XB(3)*XB(3)+RO)
+       DO 1 N=1,3
+       X(N)=XB(N)/RV
+       GRAD(N)=Q0
+   1   CONTINUE
+       RO=DSQRT(RO)
+       IF (RO.GT.EPS) GOTO 2
+       COSL=Q1
+       SINL=Q0
+       GO TO 3
+   2   CONTINUE
+       COSL=XB(1)/RO
+       SINL=XB(2)/RO
+   3   CONTINUE
+       RO=RO/RV
+       R=RR/RV
+       U=Q0
+       UU=Q0
+       U1=Q0
+       VV=Q0
+       COSK1L=Q1
+       SINK1L=Q0
+       COSKL=COSL
+       SINKL=SINL
+       PKK=Q1
+       NPN=2
+       K=1
+       NPK1=1
+       IF (DABS(RMU).LT.EPS) GOTO 100
+       V=Q0
+       RK=R
+  10   CONTINUE
+       VK=Q0
+       UK=Q0
+       RN=RK
+       P1=Q0
+       P=PKK
+       N=K
+       N1=N+1
+       NPK=NPN
+  20   CONTINUE
+       DU=RN*P
+       U1=(C(NPK1)*COSK1L+C(NPK1+1)*SINK1L)*DU+U1
+       VK=(C(NPK)*SINKL-C(NPK+1)*COSKL)*DU+VK
+       DU=(C(NPK)*COSKL+C(NPK+1)*SINKL)*DU
+       UK=DU+UK
+       V=N1*DU+V
+       IF (N.EQ.NN) GOTO 21
+       RN=RN*R
+       DU=P1
+       P1=P
+       N21=N+N1
+       P=(N21*X(3)*P1-(N+K)*DU)/(N1-K)
+       N=N1
+       N1=N+1
+       NPK1=NPK1+N21
+       NPK=NPK+N21
+       GO TO 20
+  21   CONTINUE
+       VV=K*VK+VV
+       UU=K*UK+UU
+       U=UK+U
+       IF (K.EQ.NN) GOTO 11
+       RK=RK*R
+       COSK1L=COSKL
+       SINK1L=SINKL
+       COSKL=COSK1L*COSL-SINK1L*SINL
+       SINKL=COSK1L*SINL+SINK1L*COSL
+       N21=K+K+1
+       PKK=N21*RO*PKK
+       NPK1=N21+NPN
+       NPN=NPK1+2
+       K=K+1
+       GO TO 10
+  11   CONTINUE
+       U=U*RO
+       V=V*RO
+       RN=R
+       P1=Q1
+       P=X(3)
+       N1=2
+       N=1
+       NPN=1
+  50   CONTINUE
+       DU=C(NPN)*P*RN
+       U=DU+U
+       V=N1*DU+V
+       IF (N.EQ.NN) GOTO 51
+       RN=RN*R
+       DU=P1
+       P1=P
+       N21=N+N1
+       P=(N21*X(3)*P1-N*DU)/N1
+       N=N1
+       N1=N+1
+       NPN=NPN+N21
+       GO TO 50
+  51   CONTINUE
+       DU=RMU/RV
+       U=DU*U
+       DU=-DU/RV
+       UU=U1*RO-UU*X(3)
+       GRAD(3)=(V*X(3)-UU*RO)*DU
+       UU=UU*X(3)
+       GRAD(1)=(UU*COSL-VV*SINL+V*X(1))*DU
+       GRAD(2)=(UU*SINL+VV*COSL+V*X(2))*DU
+       RETURN
+  100  CONTINUE
+       VK=Q0
+       UK=Q0
+       P1=Q0
+       P=PKK
+       N=K
+       N1=N+1
+       NPK=NPN
+  120  CONTINUE
+       U1=(C(NPK1)*COSK1L+C(NPK1+1)*SINK1L)*P+U1
+       VK=(C(NPK)*SINKL-C(NPK+1)*COSKL)*P+VK
+       UK=(C(NPK)*COSKL+C(NPK+1)*SINKL)*P+UK
+       IF (N.EQ.NN) GOTO 121
+       DU=P1
+       P1=P
+       N21=N+N1
+       P=(N21*X(3)*P1-(N+K)*DU)/(N1-K)
+       N=N1
+       N1=N+1
+       NPK1=NPK1+N21
+       NPK=NPK+N21
+       GO TO 120
+  121  CONTINUE
+       VV=K*VK+VV
+       UU=K*UK+UU
+       U=UK+U
+       IF (K.EQ.NN) GOTO 111
+       COSK1L=COSKL
+       SINK1L=SINKL
+       COSKL=COSK1L*COSL-SINK1L*SINL
+       SINKL=COSK1L*SINL+SINK1L*COSL
+       N21=K+K+1
+       PKK=N21*RO*PKK
+       NPK1=N21+NPN
+       NPN=NPK1+2
+       K=K+1
+       GO TO 100
+  111  CONTINUE
+       U=U*RO
+       P1=Q1
+       P=X(3)
+       N1=2
+       N=1
+       NPN=1
+  150  CONTINUE
+       U=C(NPN)*P+U
+       IF (N.EQ.NN) GOTO 151
+       DU=P1
+       P1=P
+       N21=N+N1
+       P=(N21*X(3)*P1-N*DU)/N1
+       N=N1
+       N1=N+1
+       NPN=NPN+N21
+       GO TO 150
+  151  CONTINUE
+       U=U*RR
+       UU=UU*X(3)-U1*RO
+       GRAD(3)=-UU*RO*R
+       UU=UU*X(3)
+       GRAD(1)=(UU*COSL+VV*SINL)*R
+       GRAD(2)=(UU*SINL-VV*COSL)*R
+       RETURN
+       END

zve/mod/nsquar.f95

+! ---------------------------------------------------------------
+!  Вычисление площади пересечения двух прямоугольников
+!  на развертке области ГСО 360x40
+!  Прямоугольник  на плоскости описывается массивом из 4 элементов
+!  A(1) - начальная долгота
+!  A(2) - конечная долгота причем A(2) >= A(1)
+!  A(3) - начальная широта
+!  A(4) - конечная широта причем A(4) >= A(3)
+!  Приведенная долгота прямоугольника  может принимать значения 0-540 град (0;360)
+!  Приведенная широта прямоугольника может принимать значения от 0 до  40 (-20;+20)
+! ---------------------------------------------------------------
+!  Bxoдныe пapaмeтpы: 
+!  Функция возвращает площадь пересечения двух прямоугольников
+!  Если возвращается 0 то пересечения нет
+!  Если чило отрицательное то вычисление невозможно
+! ---------------------------------------------------------------
+REAL*8 FUNCTION nsquar(A,B)
+    IMPLICIT NONE
+
+    REAL*8  A(4), B(4)             ! описание прямоугольников
+    REAL*8  A1,A2,A3,A4            ! Короткие обозначения сторон прямоугольника А
+    REAL*8  B1,B2,B3,B4            ! Короткие обозначения сторон прямоугольника B
+    REAL*8  C,D                    ! Стороны прямоугольника пересечения
+    integer i
+
+    nsquar=-1.0D00
+    !  Проверка правильности входных данных
+    do i=1,4                                            !  Возмодные комбинации:
+        if(A(i).lt.0.0D00.or.A(i).gt.540.0D00) return   !         A1      A2
+        if(B(i).lt.0.0D00.or.B(i).gt.540.0D00) return   !         +       +
+    enddo                                               !    +  +
+    do i=3,4                                            !1)  B1 B2                      - нет пересечения
+        if(A(i).gt.40.0D00) return                      !
+        if(B(i).gt.40.0D00) return                      !2)     +    +                  - есть пересечение
+    enddo                                               !       B1   B2
+    if(A(1).gt.A(2)) return                             !
+    if(A(3).gt.A(4)) return                             !            +  +
+    if(B(1).gt.B(2)) return                             !3)          B1 B2              - есть пересечение
+    if(B(3).gt.B(4)) return                             !               +     +
+    !  Вычисляем площадь пересечения                    !4)             B1    B2        - есть пересечение
+    nsquar=0.0D00                                       !
+    !  По долготе                                       !                     +   +
+    A1=A(1)                                             !5)                   B1  B2    - нет пересечения
+    A2=A(2)                                             !
+    B1=B(1)                                             !       +             +         - есть пересечения
+    B2=B(2)                                             !6)     B1            B2
+    D=0.0D00
+    if(B2.le.A1) return
+    if(B1.ge.A2) return
+    if(A1.ge.B1.and.A1.le.B2.and.A2.ge.B2) D=B2-A1
+    if(B1.ge.A1.and.A2.ge.B1.and.B2.ge.A2) D=A2-B1
+    if(A1.ge.B1.and.B2.ge.A2) D=A2-A1
+    if(B1.ge.A1.and.A2.ge.B2) D=B2-B1
+    !  По широте
+    A3=A(3)
+    A4=A(4)
+    B3=B(3)
+    B4=B(4)
+    C=0.0D00
+    if(B4.le.A3) return
+    if(B3.ge.A4) return
+    if(A3.ge.B3.and.A3.le.B4.and.A4.ge.B4) C=B4-A3
+    if(B3.ge.A3.and.A4.ge.B3.and.B4.ge.A4) C=A4-B3
+    if(A3.ge.B3.and.B4.ge.A4) C=A4-A3
+    if(B3.ge.A3.and.A4.ge.B4) C=B4-B3
+    !  результат с проверкой
+    if(D.lt.0.0D00.or.C.lt.0.0D00) then
+        nsquar=-2.0D00
+    else
+        nsquar=C*D
+    endif
+
+    return
+end

zve/modules/.gitignore

+!Doxyfile
+write_buildinfo.f95
\ No newline at end of file

zve/modules/makefile

+MAIN=library
+include ../makefile.inc
\ No newline at end of file