MODULE newuoa_module

    !
    !   M. J. D. Powell, "The NEWUOA software for unconstrained optimization
    !   without derivatives", in Large-Scale Nonlinear Optimization, Series:
    !   Nonconvex Optimization and Its Applications, Vol. 83, Di Pillo,
    !   Gianni; Roma, Massimo (Eds.) 2006, New York: Springer US.
    !

    ! NOTE: It is very important to not import shared_containers due to a
    ! namespace conflict fot the variable TAU that is defined here. Thus
    ! we cannot simpy use shared_interface as in all other modules.

    !/* external modules  */

    USE shared_interfaces

    USE shared_constants

    !/* setup */

    PRIVATE

    PUBLIC:: newuoa

CONTAINS
!******************************************************************************
!******************************************************************************
SUBROUTINE newuoa(FUNC, X, NPT, RHOBEG, RHOEND, IPRINT, MAXFUN, SUCCESS, MESSAGE)

    !/* external objects    */

    REAL(our_dble), INTENT(INOUT)   :: X(:)

    INTEGER(our_int), INTENT(IN)    :: NPT
    INTEGER(our_int), INTENT(IN)    :: IPRINT
    INTEGER(our_int), INTENT(IN)    :: MAXFUN

    REAL(our_dble), INTENT(IN)      :: RHOBEG
    REAL(our_dble), INTENT(IN)      :: RHOEND

    LOGICAL, INTENT(OUT)            :: SUCCESS

    CHARACTER(150), INTENT(OUT)     :: MESSAGE

    PROCEDURE(interface_func)       :: FUNC

    !/* internal objects    */

    REAL(our_dble)                  :: W(20000)

    INTEGER(our_int)                :: N

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------

    N = SIZE(X)

    CALL NEWUOA_ORIGINAL (FUNC, X, N, NPT, RHOBEG, RHOEND, IPRINT, MAXFUN, SUCCESS, MESSAGE, NF, W)

END SUBROUTINE
!******************************************************************************
!******************************************************************************
SUBROUTINE NEWUOA_ORIGINAL (FUNC, X, N, NPT, RHOBEG, RHOEND, IPRINT, MAXFUN, SUCCESS, MESSAGE, NF, W)

    IMPLICIT REAL(our_dble) (A-H,O-Z)
    DIMENSION W(*)

    REAL(our_dble)                  :: X(:)
    LOGICAL                         :: SUCCESS
    CHARACTER(150)                  :: MESSAGE

    PROCEDURE(interface_func)       :: FUNC

    !/*  internal objects       */

    INTEGER                         :: N

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     This subroutine seeks the least value of a function of many variables,
!C     by a trust region method that forms quadratic models by interpolation.
!C     There can be some freedom in the interpolation conditions, which is
!C     taken up by minimizing the Frobenius norm of the change to the second
!C     derivative of the quadratic model, beginning with a zero matrix. The
!C     arguments of the subroutine are as follows.
!C
!C     N must be set to the number of variables and must be at least two.
!C     NPT is the number of interpolation conditions. Its value must be in the
!C       interval [N+2,(N+1)(N+2)/2].
!C     Initial values of the variables must be set in X(1),X(2),...,X(N). They
!C       will be changed to the values that give the least calculated F.
!C     RHOBEG and RHOEND must be set to the initial and final values of a trust
!C       region radius, so both must be positive with RHOEND<=RHOBEG. Typically
!C       RHOBEG should be about one tenth of the greatest expected change to a
!C       variable, and RHOEND should indicate the accuracy that is required in
!C       the final values of the variables.
!C     The value of IPRINT should be set to 0, 1, 2 or 3, which controls the
!C       amount of printing. Specifically, there is no output if IPRINT=0 and
!C       there is output only at the return if IPRINT=1. Otherwise, each new
!C       value of RHO is printed, with the best vector of variables so far and
!C       the corresponding value of the objective function. Further, each new
!C       value of F with its variables are output if IPRINT=3.
!C     MAXFUN must be set to an upper bound on the number of calls of CALFUN.
!C     The array W will be used for working space. Its length must be at least
!C     (NPT+13)*(NPT+N)+3*N*(N+3)/2.
!C
!C     SUBROUTINE CALFUN (N,X,F) must be provided by the user. It must set F to
!C     the value of the objective function for the variables X(1),X(2),...,X(N).
!C
!C     Partition the working space array, so that different parts of it can be
!C     treated separately by the subroutine that performs the main calculation.
!C

    MESSAGE = 'Successful return from NEWUOA.'
    SUCCESS = .True.

    NP=N+1
    NPTM=NPT-NP
    IF (NPT .LT. N+2 .OR. NPT .GT. ((N+2)*NP)/2) THEN
      MESSAGE = 'Return from NEWUOA because NPT is not in the required interval.'
      SUCCESS = .False.
      GO TO 20
    END IF
    NDIM=NPT+N
    IXB=1
    IXO=IXB+N
    IXN=IXO+N
    IXP=IXN+N
    IFV=IXP+N*NPT
    IGQ=IFV+NPT
    IHQ=IGQ+N
    IPQ=IHQ+(N*NP)/2
    IBMAT=IPQ+NPT
    IZMAT=IBMAT+NDIM*N
    ID=IZMAT+NPT*NPTM
    IVL=ID+N
    IW=IVL+NDIM

!C
!C     The above settings provide a partition of W for subroutine NEWUOB.
!C     The partition requires the first NPT*(NPT+N)+5*N*(N+3)/2 elements of
!C     W plus the space that is needed by the last array of NEWUOB.
!C
    CALL NEWUOB (NF,FUNC,N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,W(IXB),W(IXO),W(IXN),W(IXP),W(IFV),W(IGQ),W(IHQ),W(IPQ),W(IBMAT), W(IZMAT),NDIM,W(ID),W(IVL),W(IW),SUCCESS,MESSAGE)
    20 RETURN
    END
!******************************************************************************
!******************************************************************************
    SUBROUTINE NEWUOB (NF,FUNC,N,NPT,X,RHOBEG,RHOEND,IPRINT,MAXFUN,XBASE,XOPT,XNEW,XPT,FVAL,GQ,HQ,PQ,BMAT,ZMAT,NDIM,D,VLAG,W,SUCCESS,MESSAGE)
      IMPLICIT REAL(our_dble) (A-H,O-Z)
      DIMENSION XBASE(*),XOPT(*),XNEW(*),XPT(NPT,*),FVAL(*),GQ(*),HQ(*),PQ(*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),VLAG(*),W(*)

      REAL(our_dble), INTENT(INOUT):: X(N)

      CHARACTER(150)    :: MESSAGE
      LOGICAL           :: SUCCESS

      PROCEDURE(interface_func)       :: FUNC

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     The arguments N, NPT, X, RHOBEG, RHOEND, IPRINT and MAXFUN are identical
!C       to the corresponding arguments in SUBROUTINE NEWUOA.
!C     XBASE will hold a shift of origin that should reduce the contributions
!C       from rounding errors to values of the model and Lagrange functions.
!C     XOPT will be set to the displacement from XBASE of the vector of
!C       variables that provides the least calculated F so far.
!C     XNEW will be set to the displacement from XBASE of the vector of
!C       variables for the current calculation of F.
!C     XPT will contain the interpolation point coordinates relative to XBASE.
!C     FVAL will hold the values of F at the interpolation points.
!C     GQ will hold the gradient of the quadratic model at XBASE.
!C     HQ will hold the explicit second derivatives of the quadratic model.
!C     PQ will contain the parameters of the implicit second derivatives of
!C       the quadratic model.
!C     BMAT will hold the last N columns of H.
!C     ZMAT will hold the factorization of the leading NPT by NPT submatrix of
!C       H, this factorization being ZMAT times Diag(DZ) times ZMAT^T, where
!C       the elements of DZ are plus or minus one, as specified by IDZ.
!C     NDIM is the first dimension of BMAT and has the value NPT+N.
!C     D is reserved for trial steps from XOPT.
!C     VLAG will contain the values of the Lagrange functions at a new point X.
!C       They are part of a product that requires VLAG to be of length NDIM.
!C     The array W will be used for working space. Its length must be at least
!C       10*NDIM = 10*(NPT+N).
!C
!C     Set some constants.
!C
      HALF=0.5D0
      ONE=1.0D0
      TENTH=0.1D0
      ZERO=0.0D0
      NP=N+1
      NH=(N*NP)/2
      NPTM=NPT-NP
      NFTEST=MAX0(MAXFUN,1)
!C
!C     Set the initial elements of XPT, BMAT, HQ, PQ and ZMAT to zero.
!C
      DO 20 J=1,N
      XBASE(J)=X(J)
      DO 10 K=1,NPT
   10 XPT(K,J)=ZERO
      DO 20 I=1,NDIM
   20 BMAT(I,J)=ZERO
      DO 30 IH=1,NH
   30 HQ(IH)=ZERO
      DO 40 K=1,NPT
      PQ(K)=ZERO
      DO 40 J=1,NPTM
   40 ZMAT(K,J)=ZERO
!C
!C     Begin the initialization procedure. NF becomes one more than the number
!C     of function values so far. The coordinates of the displacement of the
!C     next initial interpolation point from XBASE are set in XPT(NF,.).
!C
      RHOSQ=RHOBEG*RHOBEG
      RECIP=ONE/RHOSQ
      RECIQ=DSQRT(HALF)/RHOSQ
      NF=0
   50 NFM=NF
      NFMM=NF-N
      NF=NF+1
      IF (NFM .LE. 2*N) THEN
          IF (NFM .GE. 1 .AND. NFM .LE. N) THEN
              XPT(NF,NFM)=RHOBEG
          ELSE IF (NFM .GT. N) THEN
              XPT(NF,NFMM)=-RHOBEG
          END IF
      ELSE
          ITEMP=(NFMM-1)/N
          JPT=NFM-ITEMP*N-N
          IPT=JPT+ITEMP
          IF (IPT .GT. N) THEN
              ITEMP=JPT
              JPT=IPT-N
              IPT=ITEMP
          END IF
          XIPT=RHOBEG
          IF (FVAL(IPT+NP) .LT. FVAL(IPT+1)) XIPT=-XIPT
          XJPT=RHOBEG
          IF (FVAL(JPT+NP) .LT. FVAL(JPT+1)) XJPT=-XJPT
          XPT(NF,IPT)=XIPT
          XPT(NF,JPT)=XJPT
      END IF
!C
!C     Calculate the next value of F, label 70 being reached immediately
!C     after this calculation. The least function value so far and its index
!C     are required.
!C
      DO 60 J=1,N
   60 X(J)=XPT(NF,J)+XBASE(J)
      GOTO 310
   70 FVAL(NF)=F
      IF (NF .EQ. 1) THEN
          FBEG=F
          FOPT=F
          KOPT=1
      ELSE IF (F .LT. FOPT) THEN
          FOPT=F
          KOPT=NF
      END IF
!C
!C     Set the nonzero initial elements of BMAT and the quadratic model in
!C     the cases when NF is at most 2*N+1.
!C
      IF (NFM .LE. 2*N) THEN
          IF (NFM .GE. 1 .AND. NFM .LE. N) THEN
              GQ(NFM)=(F-FBEG)/RHOBEG
              IF (NPT .LT. NF+N) THEN
                  BMAT(1,NFM)=-ONE/RHOBEG
                  BMAT(NF,NFM)=ONE/RHOBEG
                  BMAT(NPT+NFM,NFM)=-HALF*RHOSQ
              END IF
          ELSE IF (NFM .GT. N) THEN
              BMAT(NF-N,NFMM)=HALF/RHOBEG
              BMAT(NF,NFMM)=-HALF/RHOBEG
              ZMAT(1,NFMM)=-RECIQ-RECIQ
              ZMAT(NF-N,NFMM)=RECIQ
              ZMAT(NF,NFMM)=RECIQ
              IH=(NFMM*(NFMM+1))/2
              TEMP=(FBEG-F)/RHOBEG
              HQ(IH)=(GQ(NFMM)-TEMP)/RHOBEG
              GQ(NFMM)=HALF*(GQ(NFMM)+TEMP)
          END IF
!C
!C     Set the off-diagonal second derivatives of the Lagrange functions and
!C     the initial quadratic model.
!C
      ELSE
          IH=(IPT*(IPT-1))/2+JPT
          IF (XIPT .LT. ZERO) IPT=IPT+N
          IF (XJPT .LT. ZERO) JPT=JPT+N
          ZMAT(1,NFMM)=RECIP
          ZMAT(NF,NFMM)=RECIP
          ZMAT(IPT+1,NFMM)=-RECIP
          ZMAT(JPT+1,NFMM)=-RECIP
          HQ(IH)=(FBEG-FVAL(IPT+1)-FVAL(JPT+1)+F)/(XIPT*XJPT)
      END IF
      IF (NF .LT. NPT) GOTO 50
!C
!C     Begin the iterative procedure, because the initial model is complete.
!C
      RHO=RHOBEG
      DELTA=RHO
      IDZ=1
      DIFFA=ZERO
      DIFFB=ZERO
      ITEST=0
      XOPTSQ=ZERO
      DO 80 I=1,N
      XOPT(I)=XPT(KOPT,I)
   80 XOPTSQ=XOPTSQ+XOPT(I)**2
   90 NFSAV=NF
!C
!C     Generate the next trust region step and test its length. Set KNEW
!C     to -1 if the purpose of the next F will be to improve the model.
!C
  100 KNEW=0
      CALL TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,D,W,W(NP), W(NP+N),W(NP+2*N),CRVMIN)
      DSQ=ZERO
      DO 110 I=1,N
  110 DSQ=DSQ+D(I)**2
      DNORM=DMIN1(DELTA,DSQRT(DSQ))
      IF (DNORM .LT. HALF*RHO) THEN
          KNEW=-1
          DELTA=TENTH*DELTA
          RATIO=-1.0D0
          IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
          IF (NF .LE. NFSAV+2) GOTO 460
          TEMP=0.125D0*CRVMIN*RHO*RHO
          IF (TEMP .LE. DMAX1(DIFFA,DIFFB,DIFFC)) GOTO 460
          GOTO 490
      END IF
!C
!C     Shift XBASE if XOPT may be too far from XBASE. First make the changes
!C     to BMAT that do not depend on ZMAT.
!C
  120 IF (DSQ .LE. 1.0D-3*XOPTSQ) THEN
          TEMPQ=0.25D0*XOPTSQ
          DO 140 K=1,NPT
          SUM=ZERO
          DO 130 I=1,N
  130     SUM=SUM+XPT(K,I)*XOPT(I)
          TEMP=PQ(K)*SUM
          SUM=SUM-HALF*XOPTSQ
          W(NPT+K)=SUM
          DO 140 I=1,N
          GQ(I)=GQ(I)+TEMP*XPT(K,I)
          XPT(K,I)=XPT(K,I)-HALF*XOPT(I)
          VLAG(I)=BMAT(K,I)
          W(I)=SUM*XPT(K,I)+TEMPQ*XOPT(I)
          IP=NPT+I
          DO 140 J=1,I
  140     BMAT(IP,J)=BMAT(IP,J)+VLAG(I)*W(J)+W(I)*VLAG(J)
!C
!C     Then the revisions of BMAT that depend on ZMAT are calculated.
!C
          DO 180 K=1,NPTM
          SUMZ=ZERO
          DO 150 I=1,NPT
          SUMZ=SUMZ+ZMAT(I,K)
  150     W(I)=W(NPT+I)*ZMAT(I,K)
          DO 170 J=1,N
          SUM=TEMPQ*SUMZ*XOPT(J)
          DO 160 I=1,NPT
  160     SUM=SUM+W(I)*XPT(I,J)
          VLAG(J)=SUM
          IF (K .LT. IDZ) SUM=-SUM
          DO 170 I=1,NPT
  170     BMAT(I,J)=BMAT(I,J)+SUM*ZMAT(I,K)
          DO 180 I=1,N
          IP=I+NPT
          TEMP=VLAG(I)
          IF (K .LT. IDZ) TEMP=-TEMP
          DO 180 J=1,I
  180     BMAT(IP,J)=BMAT(IP,J)+TEMP*VLAG(J)
!C
!C     The following instructions complete the shift of XBASE, including
!C     the changes to the parameters of the quadratic model.
!C
          IH=0
          DO 200 J=1,N
          W(J)=ZERO
          DO 190 K=1,NPT
          W(J)=W(J)+PQ(K)*XPT(K,J)
  190     XPT(K,J)=XPT(K,J)-HALF*XOPT(J)
          DO 200 I=1,J
          IH=IH+1
          IF (I .LT. J) GQ(J)=GQ(J)+HQ(IH)*XOPT(I)
          GQ(I)=GQ(I)+HQ(IH)*XOPT(J)
          HQ(IH)=HQ(IH)+W(I)*XOPT(J)+XOPT(I)*W(J)
  200     BMAT(NPT+I,J)=BMAT(NPT+J,I)
          DO 210 J=1,N
          XBASE(J)=XBASE(J)+XOPT(J)
  210     XOPT(J)=ZERO
          XOPTSQ=ZERO
      END IF
!C
!C     Pick the model step if KNEW is positive. A different choice of D
!C     may be made later, if the choice of D by BIGLAG causes substantial
!C     cancellation in DENOM.
!C
      IF (KNEW .GT. 0) THEN
          CALL BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DSTEP,D,ALPHA,VLAG,VLAG(NPT+1),W,W(NP),W(NP+N))
      END IF
!C
!C     Calculate VLAG and BETA for the current choice of D. The first NPT
!C     components of W_check will be held in W.
!C
      DO 230 K=1,NPT
      SUMA=ZERO
      SUMB=ZERO
      SUM=ZERO
      DO 220 J=1,N
      SUMA=SUMA+XPT(K,J)*D(J)
      SUMB=SUMB+XPT(K,J)*XOPT(J)
  220 SUM=SUM+BMAT(K,J)*D(J)
      W(K)=SUMA*(HALF*SUMA+SUMB)
  230 VLAG(K)=SUM
      BETA=ZERO
      DO 250 K=1,NPTM
      SUM=ZERO
      DO 240 I=1,NPT
  240 SUM=SUM+ZMAT(I,K)*W(I)
      IF (K .LT. IDZ) THEN
          BETA=BETA+SUM*SUM
          SUM=-SUM
      ELSE
          BETA=BETA-SUM*SUM
      END IF
      DO 250 I=1,NPT
  250 VLAG(I)=VLAG(I)+SUM*ZMAT(I,K)
      BSUM=ZERO
      DX=ZERO
      DO 280 J=1,N
      SUM=ZERO
      DO 260 I=1,NPT
  260 SUM=SUM+W(I)*BMAT(I,J)
      BSUM=BSUM+SUM*D(J)
      JP=NPT+J
      DO 270 K=1,N
  270 SUM=SUM+BMAT(JP,K)*D(K)
      VLAG(JP)=SUM
      BSUM=BSUM+SUM*D(J)
  280 DX=DX+D(J)*XOPT(J)
      BETA=DX*DX+DSQ*(XOPTSQ+DX+DX+HALF*DSQ)+BETA-BSUM
      VLAG(KOPT)=VLAG(KOPT)+ONE
!C
!C     If KNEW is positive and if the cancellation in DENOM is unacceptable,
!C     then BIGDEN calculates an alternative model step, XNEW being used for
!C     working space.
!C
      IF (KNEW .GT. 0) THEN
          TEMP=ONE+ALPHA*BETA/VLAG(KNEW)**2
          IF (DABS(TEMP) .LE. 0.8D0) THEN
              CALL BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT,KNEW,D,W,VLAG,BETA,XNEW,W(NDIM+1),W(6*NDIM+1))
          END IF
      END IF
!C
!C     Calculate the next value of the objective function.
!C
  290 DO 300 I=1,N
      XNEW(I)=XOPT(I)+D(I)
  300 X(I)=XBASE(I)+XNEW(I)
      NF=NF+1
  310 IF (NF .GT. NFTEST) THEN
          NF=NF-1
        MESSAGE = 'Return from NEWUOA because CALFUN has been called MAXFUN times.'
        SUCCESS = .False.

          GOTO 530
      END IF
      F = func(X)
      IF (IPRINT .EQ. 3) THEN
          PRINT 330, NF,F,(X(I),I=1,N)
  330      FORMAT (/4X,'Function number',I6,'    F =',1PD18.10, '    The corresponding X is:'/(2X,5D15.6))
      END IF
      IF (NF .LE. NPT) GOTO 70
      IF (KNEW .EQ. -1) GOTO 530
!C
!C     Use the quadratic model to predict the change in F due to the step D,
!C     and set DIFF to the error of this prediction.
!C
      VQUAD=ZERO
      IH=0
      DO 340 J=1,N
      VQUAD=VQUAD+D(J)*GQ(J)
      DO 340 I=1,J
      IH=IH+1
      TEMP=D(I)*XNEW(J)+D(J)*XOPT(I)
      IF (I .EQ. J) TEMP=HALF*TEMP
  340 VQUAD=VQUAD+TEMP*HQ(IH)
      DO 350 K=1,NPT
  350 VQUAD=VQUAD+PQ(K)*W(K)
      DIFF=F-FOPT-VQUAD
      DIFFC=DIFFB
      DIFFB=DIFFA
      DIFFA=DABS(DIFF)
      IF (DNORM .GT. RHO) NFSAV=NF
!C
!C     Update FOPT and XOPT if the new F is the least value of the objective
!C     function so far. The branch when KNEW is positive occurs if D is not
!C     a trust region step.
!C
      FSAVE=FOPT
      IF (F .LT. FOPT) THEN
          FOPT=F
          XOPTSQ=ZERO
          DO 360 I=1,N
          XOPT(I)=XNEW(I)
  360     XOPTSQ=XOPTSQ+XOPT(I)**2
      END IF
      KSAVE=KNEW
      IF (KNEW .GT. 0) GOTO 410
!C
!C     Pick the next value of DELTA after a trust region step.
!C
      IF (VQUAD .GE. ZERO) THEN
            MESSAGE = 'Return from NEWUOA because a trust region step has failed to reduce Q.'
            SUCCESS = .False.
          GOTO 530
      END IF
      RATIO=(F-FSAVE)/VQUAD
      IF (RATIO .LE. TENTH) THEN
          DELTA=HALF*DNORM
      ELSE IF (RATIO .LE. 0.7D0) THEN
          DELTA=DMAX1(HALF*DELTA,DNORM)
      ELSE
          DELTA=DMAX1(HALF*DELTA,DNORM+DNORM)
      END IF
      IF (DELTA .LE. 1.5D0*RHO) DELTA=RHO
!C
!C     Set KNEW to the index of the next interpolation point to be deleted.
!C
      RHOSQ=DMAX1(TENTH*DELTA,RHO)**2
      KTEMP=0
      DETRAT=ZERO
      IF (F .GE. FSAVE) THEN
          KTEMP=KOPT
          DETRAT=ONE
      END IF
      DO 400 K=1,NPT
      HDIAG=ZERO
      DO 380 J=1,NPTM
      TEMP=ONE
      IF (J .LT. IDZ) TEMP=-ONE
  380 HDIAG=HDIAG+TEMP*ZMAT(K,J)**2
      TEMP=DABS(BETA*HDIAG+VLAG(K)**2)
      DISTSQ=ZERO
      DO 390 J=1,N
  390 DISTSQ=DISTSQ+(XPT(K,J)-XOPT(J))**2
      IF (DISTSQ .GT. RHOSQ) TEMP=TEMP*(DISTSQ/RHOSQ)**3
      IF (TEMP .GT. DETRAT .AND. K .NE. KTEMP) THEN
          DETRAT=TEMP
          KNEW=K
      END IF
  400 CONTINUE
      IF (KNEW .EQ. 0) GOTO 460
!C
!C     Update BMAT, ZMAT and IDZ, so that the KNEW-th interpolation point
!C     can be moved. Begin the updating of the quadratic model, starting
!C     with the explicit second derivative term.
!C
  410 CALL UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
      FVAL(KNEW)=F
      IH=0
      DO 420 I=1,N
      TEMP=PQ(KNEW)*XPT(KNEW,I)
      DO 420 J=1,I
      IH=IH+1
  420 HQ(IH)=HQ(IH)+TEMP*XPT(KNEW,J)
      PQ(KNEW)=ZERO
!C
!C     Update the other second derivative parameters, and then the gradient
!C     vector of the model. Also include the new interpolation point.
!C
      DO 440 J=1,NPTM
      TEMP=DIFF*ZMAT(KNEW,J)
      IF (J .LT. IDZ) TEMP=-TEMP
      DO 440 K=1,NPT
  440 PQ(K)=PQ(K)+TEMP*ZMAT(K,J)
      GQSQ=ZERO
      DO 450 I=1,N
      GQ(I)=GQ(I)+DIFF*BMAT(KNEW,I)
      GQSQ=GQSQ+GQ(I)**2
  450 XPT(KNEW,I)=XNEW(I)
!C
!C     If a trust region step makes a small change to the objective function,
!C     then calculate the gradient of the least Frobenius norm interpolant at
!C     XBASE, and store it in W, using VLAG for a vector of right hand sides.
!C
      IF (KSAVE .EQ. 0 .AND. DELTA .EQ. RHO) THEN
          IF (DABS(RATIO) .GT. 1.0D-2) THEN
              ITEST=0
          ELSE
              DO 700 K=1,NPT
  700         VLAG(K)=FVAL(K)-FVAL(KOPT)
              GISQ=ZERO
              DO 720 I=1,N
              SUM=ZERO
              DO 710 K=1,NPT
  710         SUM=SUM+BMAT(K,I)*VLAG(K)
              GISQ=GISQ+SUM*SUM
  720         W(I)=SUM
!C
!C     Test whether to replace the new quadratic model by the least Frobenius
!C     norm interpolant, making the replacement if the test is satisfied.
!C
              ITEST=ITEST+1
              IF (GQSQ .LT. 1.0D2*GISQ) ITEST=0
              IF (ITEST .GE. 3) THEN
                  DO 730 I=1,N
  730             GQ(I)=W(I)
                  DO 740 IH=1,NH
  740             HQ(IH)=ZERO
                  DO 760 J=1,NPTM
                  W(J)=ZERO
                  DO 750 K=1,NPT
  750             W(J)=W(J)+VLAG(K)*ZMAT(K,J)
  760             IF (J .LT. IDZ) W(J)=-W(J)
                  DO 770 K=1,NPT
                  PQ(K)=ZERO
                  DO 770 J=1,NPTM
  770             PQ(K)=PQ(K)+ZMAT(K,J)*W(J)
                  ITEST=0
              END IF
          END IF
      END IF
      IF (F .LT. FSAVE) KOPT=KNEW
!C
!C     If a trust region step has provided a sufficient decrease in F, then
!C     branch for another trust region calculation. The case KSAVE>0 occurs
!C     when the new function value was calculated by a model step.
!C
      IF (F .LE. FSAVE+TENTH*VQUAD) GOTO 100
      IF (KSAVE .GT. 0) GOTO 100
!C
!C     Alternatively, find out if the interpolation points are close enough
!C     to the best point so far.
!C
      KNEW=0
  460 DISTSQ=4.0D0*DELTA*DELTA
      DO 480 K=1,NPT
      SUM=ZERO
      DO 470 J=1,N
  470 SUM=SUM+(XPT(K,J)-XOPT(J))**2
      IF (SUM .GT. DISTSQ) THEN
          KNEW=K
          DISTSQ=SUM
      END IF
  480 CONTINUE
!C
!C     If KNEW is positive, then set DSTEP, and branch back for the next
!C     iteration, which will generate a "model step".
!C
      IF (KNEW .GT. 0) THEN
          DSTEP=DMAX1(DMIN1(TENTH*DSQRT(DISTSQ),HALF*DELTA),RHO)
          DSQ=DSTEP*DSTEP
          GOTO 120
      END IF
      IF (RATIO .GT. ZERO) GOTO 100
      IF (DMAX1(DELTA,DNORM) .GT. RHO) GOTO 100
!C
!C     The calculations with the current value of RHO are complete. Pick the
!C     next values of RHO and DELTA.
!C
  490 IF (RHO .GT. RHOEND) THEN
          DELTA=HALF*RHO
          RATIO=RHO/RHOEND
          IF (RATIO .LE. 16.0D0) THEN
              RHO=RHOEND
          ELSE IF (RATIO .LE. 250.0D0) THEN
              RHO=DSQRT(RATIO)*RHOEND
          ELSE
              RHO=TENTH*RHO
          END IF
          DELTA=DMAX1(DELTA,RHO)
          IF (IPRINT .GE. 2) THEN
              IF (IPRINT .GE. 3) PRINT 500
  500         FORMAT (5X)
              PRINT 510, RHO,NF
  510         FORMAT (/4X,'New RHO =',1PD11.4,5X,'Number of function values =',I6)
              PRINT 520, FOPT,(XBASE(I)+XOPT(I),I=1,N)
  520         FORMAT (4X,'Least value of F =',1PD23.15,9X, 'The corresponding X is:'/(2X,5D15.6))
          END IF
          GOTO 90
      END IF
!C
!C     Return from the calculation, after another Newton-Raphson step, if
!C     it is too short to have been tried before.
!C
      IF (KNEW .EQ. -1) GOTO 290
  530 IF (FOPT .LE. F) THEN
          DO 540 I=1,N
  540     X(I)=XBASE(I)+XOPT(I)
          F=FOPT
      END IF
      IF (IPRINT .GE. 1) THEN
          PRINT 550, NF
  550     FORMAT (/4X,'At the return from NEWUOA',5X, 'Number of function values =',I6)
          PRINT 520, F,(X(I),I=1,N)
      END IF
      RETURN
      END
!******************************************************************************
!******************************************************************************
      SUBROUTINE UPDATE (N,NPT,BMAT,ZMAT,IDZ,NDIM,VLAG,BETA,KNEW,W)
      IMPLICIT REAL(our_dble) (A-H,O-Z)
      DIMENSION BMAT(NDIM,*),ZMAT(NPT,*),VLAG(*),W(*)

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     The arrays BMAT and ZMAT with IDZ are updated, in order to shift the
!C     interpolation point that has index KNEW. On entry, VLAG contains the
!C     components of the vector Theta*Wcheck+e_b of the updating formula
!C     (6.11), and BETA holds the value of the parameter that has this name.
!C     The vector W is used for working space.
!C
!C     Set some constants.
!C
      ONE=1.0D0
      ZERO=0.0D0
      NPTM=NPT-N-1
!C
!C     Apply the rotations that put zeros in the KNEW-th row of ZMAT.
!C
      JL=1
      DO 20 J=2,NPTM
      IF (J .EQ. IDZ) THEN
          JL=IDZ
      ELSE IF (ZMAT(KNEW,J) .NE. ZERO) THEN
          TEMP=DSQRT(ZMAT(KNEW,JL)**2+ZMAT(KNEW,J)**2)
          TEMPA=ZMAT(KNEW,JL)/TEMP
          TEMPB=ZMAT(KNEW,J)/TEMP
          DO 10 I=1,NPT
          TEMP=TEMPA*ZMAT(I,JL)+TEMPB*ZMAT(I,J)
          ZMAT(I,J)=TEMPA*ZMAT(I,J)-TEMPB*ZMAT(I,JL)
   10     ZMAT(I,JL)=TEMP
          ZMAT(KNEW,J)=ZERO
      END IF
   20 CONTINUE
!C
!C     Put the first NPT components of the KNEW-th column of HLAG into W,
!C     and calculate the parameters of the updating formula.
!C
      TEMPA=ZMAT(KNEW,1)
      IF (IDZ .GE. 2) TEMPA=-TEMPA
      IF (JL .GT. 1) TEMPB=ZMAT(KNEW,JL)
      DO 30 I=1,NPT
      W(I)=TEMPA*ZMAT(I,1)
      IF (JL .GT. 1) W(I)=W(I)+TEMPB*ZMAT(I,JL)
   30 CONTINUE
      ALPHA=W(KNEW)
      TAU=VLAG(KNEW)
      TAUSQ=TAU*TAU
      DENOM=ALPHA*BETA+TAUSQ
      VLAG(KNEW)=VLAG(KNEW)-ONE
!C
!C     Complete the updating of ZMAT when there is only one nonzero element
!C     in the KNEW-th row of the new matrix ZMAT, but, if IFLAG is set to one,
!C     then the first column of ZMAT will be exchanged with another one later.
!C
      IFLAG=0
      IF (JL .EQ. 1) THEN
          TEMP=DSQRT(DABS(DENOM))
          TEMPB=TEMPA/TEMP
          TEMPA=TAU/TEMP
          DO 40 I=1,NPT
   40     ZMAT(I,1)=TEMPA*ZMAT(I,1)-TEMPB*VLAG(I)
          IF (IDZ .EQ. 1 .AND. TEMP .LT. ZERO) IDZ=2
          IF (IDZ .GE. 2 .AND. TEMP .GE. ZERO) IFLAG=1
      ELSE
!C
!C     Complete the updating of ZMAT in the alternative case.
!C
          JA=1
          IF (BETA .GE. ZERO) JA=JL
          JB=JL+1-JA
          TEMP=ZMAT(KNEW,JB)/DENOM
          TEMPA=TEMP*BETA
          TEMPB=TEMP*TAU
          TEMP=ZMAT(KNEW,JA)
          SCALA=ONE/DSQRT(DABS(BETA)*TEMP*TEMP+TAUSQ)
          SCALB=SCALA*DSQRT(DABS(DENOM))
          DO 50 I=1,NPT
          ZMAT(I,JA)=SCALA*(TAU*ZMAT(I,JA)-TEMP*VLAG(I))
   50     ZMAT(I,JB)=SCALB*(ZMAT(I,JB)-TEMPA*W(I)-TEMPB*VLAG(I))
          IF (DENOM .LE. ZERO) THEN
              IF (BETA .LT. ZERO) IDZ=IDZ+1
              IF (BETA .GE. ZERO) IFLAG=1
          END IF
      END IF
!C
!C     IDZ is reduced in the following case, and usually the first column
!C     of ZMAT is exchanged with a later one.
!C
      IF (IFLAG .EQ. 1) THEN
          IDZ=IDZ-1
          DO 60 I=1,NPT
          TEMP=ZMAT(I,1)
          ZMAT(I,1)=ZMAT(I,IDZ)
   60     ZMAT(I,IDZ)=TEMP
      END IF
!C
!C     Finally, update the matrix BMAT.
!C
      DO 70 J=1,N
      JP=NPT+J
      W(JP)=BMAT(KNEW,J)
      TEMPA=(ALPHA*VLAG(JP)-TAU*W(JP))/DENOM
      TEMPB=(-BETA*W(JP)-TAU*VLAG(JP))/DENOM
      DO 70 I=1,JP
      BMAT(I,J)=BMAT(I,J)+TEMPA*VLAG(I)+TEMPB*W(I)
      IF (I .GT. NPT) BMAT(JP,I-NPT)=BMAT(I,J)
   70 CONTINUE
      RETURN
      END
!******************************************************************************
!******************************************************************************
      SUBROUTINE TRSAPP (N,NPT,XOPT,XPT,GQ,HQ,PQ,DELTA,STEP, D,G,HD,HS,CRVMIN)
      IMPLICIT REAL(our_dble) (A-H,O-Z)
      DIMENSION XOPT(*),XPT(NPT,*),GQ(*),HQ(*),PQ(*),STEP(*), D(*),G(*),HD(*),HS(*)

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     N is the number of variables of a quadratic objective function, Q say.
!C     The arguments NPT, XOPT, XPT, GQ, HQ and PQ have their usual meanings,
!C       in order to define the current quadratic model Q.
!C     DELTA is the trust region radius, and has to be positive.
!C     STEP will be set to the calculated trial step.
!C     The arrays D, G, HD and HS will be used for working space.
!C     CRVMIN will be set to the least curvature of H along the conjugate
!C       directions that occur, except that it is set to zero if STEP goes
!C       all the way to the trust region boundary.
!C
!C     The calculation of STEP begins with the truncated conjugate gradient
!C     method. If the boundary of the trust region is reached, then further
!C     changes to STEP may be made, each one being in the 2D space spanned
!C     by the current STEP and the corresponding gradient of Q. Thus STEP
!C     should provide a substantial reduction to Q within the trust region.
!C
!C     Initialization, which includes setting HD to H times XOPT.
!C
      HALF=0.5D0
      ZERO=0.0D0
      TWOPI=8.0D0*DATAN(1.0D0)
      DELSQ=DELTA*DELTA
      ITERC=0
      ITERMAX=N
      ITERSW=ITERMAX
      DO 10 I=1,N
   10 D(I)=XOPT(I)
      GOTO 170
!C
!C     Prepare for the first line search.
!C
   20 QRED=ZERO
      DD=ZERO
      DO 30 I=1,N
      STEP(I)=ZERO
      HS(I)=ZERO
      G(I)=GQ(I)+HD(I)
      D(I)=-G(I)
   30 DD=DD+D(I)**2
      CRVMIN=ZERO
      IF (DD .EQ. ZERO) GOTO 160
      DS=ZERO
      SS=ZERO
      GG=DD
      GGBEG=GG
!C
!C     Calculate the step to the trust region boundary and the product HD.
!C
   40 ITERC=ITERC+1
      TEMP=DELSQ-SS
      BSTEP=TEMP/(DS+DSQRT(DS*DS+DD*TEMP))
      GOTO 170
   50 DHD=ZERO
      DO 60 J=1,N
   60 DHD=DHD+D(J)*HD(J)
!C
!C     Update CRVMIN and set the step-length ALPHA.
!C
      ALPHA=BSTEP
      IF (DHD .GT. ZERO) THEN
          TEMP=DHD/DD
          IF (ITERC .EQ. 1) CRVMIN=TEMP
          CRVMIN=DMIN1(CRVMIN,TEMP)
          ALPHA=DMIN1(ALPHA,GG/DHD)
      END IF
      QADD=ALPHA*(GG-HALF*ALPHA*DHD)
      QRED=QRED+QADD
!C
!C     Update STEP and HS.
!C
      GGSAV=GG
      GG=ZERO
      DO 70 I=1,N
      STEP(I)=STEP(I)+ALPHA*D(I)
      HS(I)=HS(I)+ALPHA*HD(I)
   70 GG=GG+(G(I)+HS(I))**2
!C
!C     Begin another conjugate direction iteration if required.
!C
      IF (ALPHA .LT. BSTEP) THEN
          IF (QADD .LE. 0.01D0*QRED) GOTO 160
          IF (GG .LE. 1.0D-4*GGBEG) GOTO 160
          IF (ITERC .EQ. ITERMAX) GOTO 160
          TEMP=GG/GGSAV
          DD=ZERO
          DS=ZERO
          SS=ZERO
          DO 80 I=1,N
          D(I)=TEMP*D(I)-G(I)-HS(I)
          DD=DD+D(I)**2
          DS=DS+D(I)*STEP(I)
   80     SS=SS+STEP(I)**2
          IF (DS .LE. ZERO) GOTO 160
          IF (SS .LT. DELSQ) GOTO 40
      END IF
      CRVMIN=ZERO
      ITERSW=ITERC
!C
!C     Test whether an alternative iteration is required.
!C
   90 IF (GG .LE. 1.0D-4*GGBEG) GOTO 160
      SG=ZERO
      SHS=ZERO
      DO 100 I=1,N
      SG=SG+STEP(I)*G(I)
  100 SHS=SHS+STEP(I)*HS(I)
      SGK=SG+SHS
      ANGTEST=SGK/DSQRT(GG*DELSQ)
      IF (ANGTEST .LE. -0.99D0) GOTO 160
!C
!C     Begin the alternative iteration by calculating D and HD and some
!C     scalar products.
!C
      ITERC=ITERC+1
      TEMP=DSQRT(DELSQ*GG-SGK*SGK)
      TEMPA=DELSQ/TEMP
      TEMPB=SGK/TEMP
      DO 110 I=1,N
  110 D(I)=TEMPA*(G(I)+HS(I))-TEMPB*STEP(I)
      GOTO 170
  120 DG=ZERO
      DHD=ZERO
      DHS=ZERO
      DO 130 I=1,N
      DG=DG+D(I)*G(I)
      DHD=DHD+HD(I)*D(I)
  130 DHS=DHS+HD(I)*STEP(I)
!C
!C     Seek the value of the angle that minimizes Q.
!C
      CF=HALF*(SHS-DHD)
      QBEG=SG+CF
      QSAV=QBEG
      QMIN=QBEG
      ISAVE=0
      IU=49
      TEMP=TWOPI/DBLE(IU+1)
      DO 140 I=1,IU
      ANGLE=DBLE(I)*TEMP
      CTH=DCOS(ANGLE)
      STH=DSIN(ANGLE)
      QNEW=(SG+CF*CTH)*CTH+(DG+DHS*CTH)*STH
      IF (QNEW .LT. QMIN) THEN
          QMIN=QNEW
          ISAVE=I
          TEMPA=QSAV
      ELSE IF (I .EQ. ISAVE+1) THEN
          TEMPB=QNEW
      END IF
  140 QSAV=QNEW
      IF (ISAVE .EQ. ZERO) TEMPA=QNEW
      IF (ISAVE .EQ. IU) TEMPB=QBEG
      ANGLE=ZERO
      IF (TEMPA .NE. TEMPB) THEN
          TEMPA=TEMPA-QMIN
          TEMPB=TEMPB-QMIN
          ANGLE=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
      END IF
      ANGLE=TEMP*(DBLE(ISAVE)+ANGLE)
!C
!C     Calculate the new STEP and HS. Then test for convergence.
!C
      CTH=DCOS(ANGLE)
      STH=DSIN(ANGLE)
      REDUC=QBEG-(SG+CF*CTH)*CTH-(DG+DHS*CTH)*STH
      GG=ZERO
      DO 150 I=1,N
      STEP(I)=CTH*STEP(I)+STH*D(I)
      HS(I)=CTH*HS(I)+STH*HD(I)
  150 GG=GG+(G(I)+HS(I))**2
      QRED=QRED+REDUC
      RATIO=REDUC/QRED
      IF (ITERC .LT. ITERMAX .AND. RATIO .GT. 0.01D0) GOTO 90
  160 RETURN
!C
!C     The following instructions act as a subroutine for setting the vector
!C     HD to the vector D multiplied by the second derivative matrix of Q.
!C     They are called from three different places, which are distinguished
!C     by the value of ITERC.
!C
  170 DO 180 I=1,N
  180 HD(I)=ZERO
      DO 200 K=1,NPT
      TEMP=ZERO
      DO 190 J=1,N
  190 TEMP=TEMP+XPT(K,J)*D(J)
      TEMP=TEMP*PQ(K)
      DO 200 I=1,N
  200 HD(I)=HD(I)+TEMP*XPT(K,I)
      IH=0
      DO 210 J=1,N
      DO 210 I=1,J
      IH=IH+1
      IF (I .LT. J) HD(J)=HD(J)+HQ(IH)*D(I)
  210 HD(I)=HD(I)+HQ(IH)*D(J)
      IF (ITERC .EQ. 0) GOTO 20
      IF (ITERC .LE. ITERSW) GOTO 50
      GOTO 120
      END
!******************************************************************************
!******************************************************************************
      SUBROUTINE BIGLAG (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KNEW,DELTA,D,ALPHA,HCOL,GC,GD,S,W)
      IMPLICIT REAL(our_dble) (A-H,O-Z)
      DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),HCOL(*),GC(*),GD(*),S(*),W(*)

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     N is the number of variables.
!C     NPT is the number of interpolation equations.
!C     XOPT is the best interpolation point so far.
!C     XPT contains the coordinates of the current interpolation points.
!C     BMAT provides the last N columns of H.
!C     ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H.
!C     NDIM is the first dimension of BMAT and has the value NPT+N.
!C     KNEW is the index of the interpolation point that is going to be moved.
!C     DELTA is the current trust region bound.
!C     D will be set to the step from XOPT to the new point.
!C     ALPHA will be set to the KNEW-th diagonal element of the H matrix.
!C     HCOL, GC, GD, S and W will be used for working space.
!C
!C     The step D is calculated in a way that attempts to maximize the modulus
!C     of LFUNC(XOPT+D), subject to the bound ||D|| .LE. DELTA, where LFUNC is
!C     the KNEW-th Lagrange function.
!C
!C     Set some constants.
!C
      HALF=0.5D0
      ONE=1.0D0
      ZERO=0.0D0
      TWOPI=8.0D0*DATAN(ONE)
      DELSQ=DELTA*DELTA
      NPTM=NPT-N-1
!C
!C     Set the first NPT components of HCOL to the leading elements of the
!C     KNEW-th column of H.
!C
      ITERC=0
      DO 10 K=1,NPT
   10 HCOL(K)=ZERO
      DO 20 J=1,NPTM
      TEMP=ZMAT(KNEW,J)
      IF (J .LT. IDZ) TEMP=-TEMP
      DO 20 K=1,NPT
   20 HCOL(K)=HCOL(K)+TEMP*ZMAT(K,J)
      ALPHA=HCOL(KNEW)
!C
!C     Set the unscaled initial direction D. Form the gradient of LFUNC at
!C     XOPT, and multiply D by the second derivative matrix of LFUNC.
!C
      DD=ZERO
      DO 30 I=1,N
      D(I)=XPT(KNEW,I)-XOPT(I)
      GC(I)=BMAT(KNEW,I)
      GD(I)=ZERO
   30 DD=DD+D(I)**2
      DO 50 K=1,NPT
      TEMP=ZERO
      SUM=ZERO
      DO 40 J=1,N
      TEMP=TEMP+XPT(K,J)*XOPT(J)
   40 SUM=SUM+XPT(K,J)*D(J)
      TEMP=HCOL(K)*TEMP
      SUM=HCOL(K)*SUM
      DO 50 I=1,N
      GC(I)=GC(I)+TEMP*XPT(K,I)
   50 GD(I)=GD(I)+SUM*XPT(K,I)
!C
!C     Scale D and GD, with a sign change if required. Set S to another
!C     vector in the initial two dimensional subspace.
!C
      GG=ZERO
      SP=ZERO
      DHD=ZERO
      DO 60 I=1,N
      GG=GG+GC(I)**2
      SP=SP+D(I)*GC(I)
   60 DHD=DHD+D(I)*GD(I)
      SCALE=DELTA/DSQRT(DD)
      IF (SP*DHD .LT. ZERO) SCALE=-SCALE
      TEMP=ZERO
      IF (SP*SP .GT. 0.99D0*DD*GG) TEMP=ONE
      TAU=SCALE*(DABS(SP)+HALF*SCALE*DABS(DHD))
      IF (GG*DELSQ .LT. 0.01D0*TAU*TAU) TEMP=ONE
      DO 70 I=1,N
      D(I)=SCALE*D(I)
      GD(I)=SCALE*GD(I)
   70 S(I)=GC(I)+TEMP*GD(I)
!C
!C     Begin the iteration by overwriting S with a vector that has the
!C     required length and direction, except that termination occurs if
!C     the given D and S are nearly parallel.
!C
   80 ITERC=ITERC+1
      DD=ZERO
      SP=ZERO
      SS=ZERO
      DO 90 I=1,N
      DD=DD+D(I)**2
      SP=SP+D(I)*S(I)
   90 SS=SS+S(I)**2
      TEMP=DD*SS-SP*SP
      IF (TEMP .LE. 1.0D-8*DD*SS) GOTO 160
      DENOM=DSQRT(TEMP)
      DO 100 I=1,N
      S(I)=(DD*S(I)-SP*D(I))/DENOM
  100 W(I)=ZERO
!C
!C     Calculate the coefficients of the objective function on the circle,
!C     beginning with the multiplication of S by the second derivative matrix.
!C
      DO 120 K=1,NPT
      SUM=ZERO
      DO 110 J=1,N
  110 SUM=SUM+XPT(K,J)*S(J)
      SUM=HCOL(K)*SUM
      DO 120 I=1,N
  120 W(I)=W(I)+SUM*XPT(K,I)
      CF1=ZERO
      CF2=ZERO
      CF3=ZERO
      CF4=ZERO
      CF5=ZERO
      DO 130 I=1,N
      CF1=CF1+S(I)*W(I)
      CF2=CF2+D(I)*GC(I)
      CF3=CF3+S(I)*GC(I)
      CF4=CF4+D(I)*GD(I)
  130 CF5=CF5+S(I)*GD(I)
      CF1=HALF*CF1
      CF4=HALF*CF4-CF1
!C
!C     Seek the value of the angle that maximizes the modulus of TAU.
!C
      TAUBEG=CF1+CF2+CF4
      TAUMAX=TAUBEG
      TAUOLD=TAUBEG
      ISAVE=0
      IU=49
      TEMP=TWOPI/DBLE(IU+1)
      DO 140 I=1,IU
      ANGLE=DBLE(I)*TEMP
      CTH=DCOS(ANGLE)
      STH=DSIN(ANGLE)
      TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
      IF (DABS(TAU) .GT. DABS(TAUMAX)) THEN
          TAUMAX=TAU
          ISAVE=I
          TEMPA=TAUOLD
      ELSE IF (I .EQ. ISAVE+1) THEN
          TEMPB=TAU
      END IF
  140 TAUOLD=TAU
      IF (ISAVE .EQ. 0) TEMPA=TAU
      IF (ISAVE .EQ. IU) TEMPB=TAUBEG
      STEP=ZERO
      IF (TEMPA .NE. TEMPB) THEN
          TEMPA=TEMPA-TAUMAX
          TEMPB=TEMPB-TAUMAX
          STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
      END IF
      ANGLE=TEMP*(DBLE(ISAVE)+STEP)
!C
!C     Calculate the new D and GD. Then test for convergence.
!C
      CTH=DCOS(ANGLE)
      STH=DSIN(ANGLE)
      TAU=CF1+(CF2+CF4*CTH)*CTH+(CF3+CF5*CTH)*STH
      DO 150 I=1,N
      D(I)=CTH*D(I)+STH*S(I)
      GD(I)=CTH*GD(I)+STH*W(I)
  150 S(I)=GC(I)+GD(I)
      IF (DABS(TAU) .LE. 1.1D0*DABS(TAUBEG)) GOTO 160
      IF (ITERC .LT. N) GOTO 80
  160 RETURN
      END
!******************************************************************************
!******************************************************************************
      SUBROUTINE BIGDEN (N,NPT,XOPT,XPT,BMAT,ZMAT,IDZ,NDIM,KOPT,KNEW,D,W,VLAG,BETA,S,WVEC,PROD)
      IMPLICIT REAL(our_dble) (A-H,O-Z)
      DIMENSION XOPT(*),XPT(NPT,*),BMAT(NDIM,*),ZMAT(NPT,*),D(*),W(*),VLAG(*),S(*),WVEC(NDIM,*),PROD(NDIM,*)
      DIMENSION DEN(9),DENEX(9),PAR(9)

!------------------------------------------------------------------------------
! Algorithm
!------------------------------------------------------------------------------
!C
!C     N is the number of variables.
!C     NPT is the number of interpolation equations.
!C     XOPT is the best interpolation point so far.
!C     XPT contains the coordinates of the current interpolation points.
!C     BMAT provides the last N columns of H.
!C     ZMAT and IDZ give a factorization of the first NPT by NPT submatrix of H.
!C     NDIM is the first dimension of BMAT and has the value NPT+N.
!C     KOPT is the index of the optimal interpolation point.
!C     KNEW is the index of the interpolation point that is going to be moved.
!C     D will be set to the step from XOPT to the new point, and on entry it
!C       should be the D that was calculated by the last call of BIGLAG. The
!C       length of the initial D provides a trust region bound on the final D.
!C     W will be set to Wcheck for the final choice of D.
!C     VLAG will be set to Theta*Wcheck+e_b for the final choice of D.
!C     BETA will be set to the value that will occur in the updating formula
!C       when the KNEW-th interpolation point is moved to its new position.
!C     S, WVEC, PROD and the private arrays DEN, DENEX and PAR will be used
!C       for working space.
!C
!C     D is calculated in a way that should provide a denominator with a large
!C     modulus in the updating formula when the KNEW-th interpolation point is
!C     shifted to the new position XOPT+D.
!C
!C     Set some constants.
!C
      HALF=0.5D0
      ONE=1.0D0
      QUART=0.25D0
      TWO=2.0D0
      ZERO=0.0D0
      TWOPI=8.0D0*DATAN(ONE)
      NPTM=NPT-N-1
!C
!C     Store the first NPT elements of the KNEW-th column of H in W(N+1)
!C     to W(N+NPT).
!C
      DO 10 K=1,NPT
   10 W(N+K)=ZERO
      DO 20 J=1,NPTM
      TEMP=ZMAT(KNEW,J)
      IF (J .LT. IDZ) TEMP=-TEMP
      DO 20 K=1,NPT
   20 W(N+K)=W(N+K)+TEMP*ZMAT(K,J)
      ALPHA=W(N+KNEW)
!C
!C     The initial search direction D is taken from the last call of BIGLAG,
!C     and the initial S is set below, usually to the direction from X_OPT
!C     to X_KNEW, but a different direction to an interpolation point may
!C     be chosen, in order to prevent S from being nearly parallel to D.
!C
      DD=ZERO
      DS=ZERO
      SS=ZERO
      XOPTSQ=ZERO
      DO 30 I=1,N
      DD=DD+D(I)**2
      S(I)=XPT(KNEW,I)-XOPT(I)
      DS=DS+D(I)*S(I)
      SS=SS+S(I)**2
   30 XOPTSQ=XOPTSQ+XOPT(I)**2
      IF (DS*DS .GT. 0.99D0*DD*SS) THEN
          KSAV=KNEW
          DTEST=DS*DS/SS
          DO 50 K=1,NPT
          IF (K .NE. KOPT) THEN
              DSTEMP=ZERO
              SSTEMP=ZERO
              DO 40 I=1,N
              DIFF=XPT(K,I)-XOPT(I)
              DSTEMP=DSTEMP+D(I)*DIFF
   40         SSTEMP=SSTEMP+DIFF*DIFF
              IF (DSTEMP*DSTEMP/SSTEMP .LT. DTEST) THEN
                  KSAV=K
                  DTEST=DSTEMP*DSTEMP/SSTEMP
                  DS=DSTEMP
                  SS=SSTEMP
              END IF
          END IF
   50     CONTINUE
          DO 60 I=1,N
   60     S(I)=XPT(KSAV,I)-XOPT(I)
      END IF
      SSDEN=DD*SS-DS*DS
      ITERC=0
      DENSAV=ZERO
!C
!C     Begin the iteration by overwriting S with a vector that has the
!C     required length and direction.
!C
   70 ITERC=ITERC+1
      TEMP=ONE/DSQRT(SSDEN)
      XOPTD=ZERO
      XOPTS=ZERO
      DO 80 I=1,N
      S(I)=TEMP*(DD*S(I)-DS*D(I))
      XOPTD=XOPTD+XOPT(I)*D(I)
   80 XOPTS=XOPTS+XOPT(I)*S(I)
!C
!C     Set the coefficients of the first two terms of BETA.
!C
      TEMPA=HALF*XOPTD*XOPTD
      TEMPB=HALF*XOPTS*XOPTS
      DEN(1)=DD*(XOPTSQ+HALF*DD)+TEMPA+TEMPB
      DEN(2)=TWO*XOPTD*DD
      DEN(3)=TWO*XOPTS*DD
      DEN(4)=TEMPA-TEMPB
      DEN(5)=XOPTD*XOPTS
      DO 90 I=6,9
   90 DEN(I)=ZERO
!C
!C     Put the coefficients of Wcheck in WVEC.
!C
      DO 110 K=1,NPT
      TEMPA=ZERO
      TEMPB=ZERO
      TEMPC=ZERO
      DO 100 I=1,N
      TEMPA=TEMPA+XPT(K,I)*D(I)
      TEMPB=TEMPB+XPT(K,I)*S(I)
  100 TEMPC=TEMPC+XPT(K,I)*XOPT(I)
      WVEC(K,1)=QUART*(TEMPA*TEMPA+TEMPB*TEMPB)
      WVEC(K,2)=TEMPA*TEMPC
      WVEC(K,3)=TEMPB*TEMPC
      WVEC(K,4)=QUART*(TEMPA*TEMPA-TEMPB*TEMPB)
  110 WVEC(K,5)=HALF*TEMPA*TEMPB
      DO 120 I=1,N
      IP=I+NPT
      WVEC(IP,1)=ZERO
      WVEC(IP,2)=D(I)
      WVEC(IP,3)=S(I)
      WVEC(IP,4)=ZERO
  120 WVEC(IP,5)=ZERO
!C
!C     Put the coefficents of THETA*Wcheck in PROD.
!C
      DO 190 JC=1,5
      NW=NPT
      IF (JC .EQ. 2 .OR. JC .EQ. 3) NW=NDIM
      DO 130 K=1,NPT
  130 PROD(K,JC)=ZERO
      DO 150 J=1,NPTM
      SUM=ZERO
      DO 140 K=1,NPT
  140 SUM=SUM+ZMAT(K,J)*WVEC(K,JC)
      IF (J .LT. IDZ) SUM=-SUM
      DO 150 K=1,NPT
  150 PROD(K,JC)=PROD(K,JC)+SUM*ZMAT(K,J)
      IF (NW .EQ. NDIM) THEN
          DO 170 K=1,NPT
          SUM=ZERO
          DO 160 J=1,N
  160     SUM=SUM+BMAT(K,J)*WVEC(NPT+J,JC)
  170     PROD(K,JC)=PROD(K,JC)+SUM
      END IF
      DO 190 J=1,N
      SUM=ZERO
      DO 180 I=1,NW
  180 SUM=SUM+BMAT(I,J)*WVEC(I,JC)
  190 PROD(NPT+J,JC)=SUM
!C
!C     Include in DEN the part of BETA that depends on THETA.
!C
      DO 210 K=1,NDIM
      SUM=ZERO
      DO 200 I=1,5
      PAR(I)=HALF*PROD(K,I)*WVEC(K,I)
  200 SUM=SUM+PAR(I)
      DEN(1)=DEN(1)-PAR(1)-SUM
      TEMPA=PROD(K,1)*WVEC(K,2)+PROD(K,2)*WVEC(K,1)
      TEMPB=PROD(K,2)*WVEC(K,4)+PROD(K,4)*WVEC(K,2)
      TEMPC=PROD(K,3)*WVEC(K,5)+PROD(K,5)*WVEC(K,3)
      DEN(2)=DEN(2)-TEMPA-HALF*(TEMPB+TEMPC)
      DEN(6)=DEN(6)-HALF*(TEMPB-TEMPC)
      TEMPA=PROD(K,1)*WVEC(K,3)+PROD(K,3)*WVEC(K,1)
      TEMPB=PROD(K,2)*WVEC(K,5)+PROD(K,5)*WVEC(K,2)
      TEMPC=PROD(K,3)*WVEC(K,4)+PROD(K,4)*WVEC(K,3)
      DEN(3)=DEN(3)-TEMPA-HALF*(TEMPB-TEMPC)
      DEN(7)=DEN(7)-HALF*(TEMPB+TEMPC)
      TEMPA=PROD(K,1)*WVEC(K,4)+PROD(K,4)*WVEC(K,1)
      DEN(4)=DEN(4)-TEMPA-PAR(2)+PAR(3)
      TEMPA=PROD(K,1)*WVEC(K,5)+PROD(K,5)*WVEC(K,1)
      TEMPB=PROD(K,2)*WVEC(K,3)+PROD(K,3)*WVEC(K,2)
      DEN(5)=DEN(5)-TEMPA-HALF*TEMPB
      DEN(8)=DEN(8)-PAR(4)+PAR(5)
      TEMPA=PROD(K,4)*WVEC(K,5)+PROD(K,5)*WVEC(K,4)
  210 DEN(9)=DEN(9)-HALF*TEMPA
!C
!C     Extend DEN so that it holds all the coefficients of DENOM.
!C
      SUM=ZERO
      DO 220 I=1,5
      PAR(I)=HALF*PROD(KNEW,I)**2
  220 SUM=SUM+PAR(I)
      DENEX(1)=ALPHA*DEN(1)+PAR(1)+SUM
      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,2)
      TEMPB=PROD(KNEW,2)*PROD(KNEW,4)
      TEMPC=PROD(KNEW,3)*PROD(KNEW,5)
      DENEX(2)=ALPHA*DEN(2)+TEMPA+TEMPB+TEMPC
      DENEX(6)=ALPHA*DEN(6)+TEMPB-TEMPC
      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,3)
      TEMPB=PROD(KNEW,2)*PROD(KNEW,5)
      TEMPC=PROD(KNEW,3)*PROD(KNEW,4)
      DENEX(3)=ALPHA*DEN(3)+TEMPA+TEMPB-TEMPC
      DENEX(7)=ALPHA*DEN(7)+TEMPB+TEMPC
      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,4)
      DENEX(4)=ALPHA*DEN(4)+TEMPA+PAR(2)-PAR(3)
      TEMPA=TWO*PROD(KNEW,1)*PROD(KNEW,5)
      DENEX(5)=ALPHA*DEN(5)+TEMPA+PROD(KNEW,2)*PROD(KNEW,3)
      DENEX(8)=ALPHA*DEN(8)+PAR(4)-PAR(5)
      DENEX(9)=ALPHA*DEN(9)+PROD(KNEW,4)*PROD(KNEW,5)
!C
!C     Seek the value of the angle that maximizes the modulus of DENOM.
!C
      SUM=DENEX(1)+DENEX(2)+DENEX(4)+DENEX(6)+DENEX(8)
      DENOLD=SUM
      DENMAX=SUM
      ISAVE=0
      IU=49
      TEMP=TWOPI/DBLE(IU+1)
      PAR(1)=ONE
      DO 250 I=1,IU
      ANGLE=DBLE(I)*TEMP
      PAR(2)=DCOS(ANGLE)
      PAR(3)=DSIN(ANGLE)
      DO 230 J=4,8,2
      PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
  230 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
      SUMOLD=SUM
      SUM=ZERO
      DO 240 J=1,9
  240 SUM=SUM+DENEX(J)*PAR(J)
      IF (DABS(SUM) .GT. DABS(DENMAX)) THEN
          DENMAX=SUM
          ISAVE=I
          TEMPA=SUMOLD
      ELSE IF (I .EQ. ISAVE+1) THEN
          TEMPB=SUM
      END IF
  250 CONTINUE
      IF (ISAVE .EQ. 0) TEMPA=SUM
      IF (ISAVE .EQ. IU) TEMPB=DENOLD
      STEP=ZERO
      IF (TEMPA .NE. TEMPB) THEN
          TEMPA=TEMPA-DENMAX
          TEMPB=TEMPB-DENMAX
          STEP=HALF*(TEMPA-TEMPB)/(TEMPA+TEMPB)
      END IF
      ANGLE=TEMP*(DBLE(ISAVE)+STEP)
!C
!C     Calculate the new parameters of the denominator, the new VLAG vector
!C     and the new D. Then test for convergence.
!C
      PAR(2)=DCOS(ANGLE)
      PAR(3)=DSIN(ANGLE)
      DO 260 J=4,8,2
      PAR(J)=PAR(2)*PAR(J-2)-PAR(3)*PAR(J-1)
  260 PAR(J+1)=PAR(2)*PAR(J-1)+PAR(3)*PAR(J-2)
      BETA=ZERO
      DENMAX=ZERO
      DO 270 J=1,9
      BETA=BETA+DEN(J)*PAR(J)
  270 DENMAX=DENMAX+DENEX(J)*PAR(J)
      DO 280 K=1,NDIM
      VLAG(K)=ZERO
      DO 280 J=1,5
  280 VLAG(K)=VLAG(K)+PROD(K,J)*PAR(J)
      TAU=VLAG(KNEW)
      DD=ZERO
      TEMPA=ZERO
      TEMPB=ZERO
      DO 290 I=1,N
      D(I)=PAR(2)*D(I)+PAR(3)*S(I)
      W(I)=XOPT(I)+D(I)
      DD=DD+D(I)**2
      TEMPA=TEMPA+D(I)*W(I)
  290 TEMPB=TEMPB+W(I)*W(I)
      IF (ITERC .GE. N) GOTO 340
      IF (ITERC .GT. 1) DENSAV=DMAX1(DENSAV,DENOLD)
      IF (DABS(DENMAX) .LE. 1.1D0*DABS(DENSAV)) GOTO 340
      DENSAV=DENMAX
!C
!C     Set S to half the gradient of the denominator with respect to D.
!C     Then branch for the next iteration.
!C
      DO 300 I=1,N
      TEMP=TEMPA*XOPT(I)+TEMPB*D(I)-VLAG(NPT+I)
  300 S(I)=TAU*BMAT(KNEW,I)+ALPHA*TEMP
      DO 320 K=1,NPT
      SUM=ZERO
      DO 310 J=1,N
  310 SUM=SUM+XPT(K,J)*W(J)
      TEMP=(TAU*W(N+K)-ALPHA*VLAG(K))*SUM
      DO 320 I=1,N
  320 S(I)=S(I)+TEMP*XPT(K,I)
      SS=ZERO
      DS=ZERO
      DO 330 I=1,N
      SS=SS+S(I)**2
  330 DS=DS+D(I)*S(I)
      SSDEN=DD*SS-DS*DS
      IF (SSDEN .GE. 1.0D-8*DD*SS) GOTO 70
!C
!C     Set the vector W before the RETURN from the subroutine.
!C
  340 DO 350 K=1,NDIM
      W(K)=ZERO
      DO 350 J=1,5
  350 W(K)=W(K)+WVEC(K,J)*PAR(J)
      VLAG(KOPT)=VLAG(KOPT)+ONE
      RETURN
      END
!******************************************************************************
!******************************************************************************
END MODULE


text:
C
C       DUMMY DEFINITIONS FOR NCAR GRAPHICS FUNCTIONS
C
        SUBROUTINE EZXY
        END
        SUBROUTINE ANOTAT
        END
        SUBROUTINE AGSETI
        END
        SUBROUTINE GCLKS
        END
        SUBROUTINE GCLWK
        END
        SUBROUTINE GDAWK
        END
        SUBROUTINE GACWK
        END
        SUBROUTINE GOPWK
        END
        SUBROUTINE GOPKS
        END
        SUBROUTINE MAPDRW
        END
        SUBROUTINE MAPINT
        END
        SUBROUTINE MAPSET
        END
        SUBROUTINE MAPROJ
        END
        SUBROUTINE MAPSTI
        END
        SUBROUTINE MAPSTC
        END
        SUBROUTINE FRAME
        END
        SUBROUTINE CPSETC
        END
        SUBROUTINE CPSETR
        END
        SUBROUTINE CPSETI
        END
        SUBROUTINE WTSTR
        END
        SUBROUTINE Q8QST4
        END
        SUBROUTINE R1MACH
        END
        SUBROUTINE CPUY
        END
        SUBROUTINE CPUX
        END
        SUBROUTINE ENCD
        END
        SUBROUTINE GQNT
        END
        SUBROUTINE SET
        END
        SUBROUTINE DRWVEC
        END
        SUBROUTINE PERIM
        END
        SUBROUTINE FL2INT
        END
        SUBROUTINE GSELNT
        END
        SUBROUTINE GSCLIP
        END
        SUBROUTINE GQCLIP
        END
        SUBROUTINE GQCNTN
        END
        SUBROUTINE GETUSV
        END
        SUBROUTINE IOR
        END
        SUBROUTINE IAND
        END
        SUBROUTINE GSASF
        END
        SUBROUTINE ARSCAM
        END
        SUBROUTINE ARINAM
        END
        SUBROUTINE CPCLAM
        END
        SUBROUTINE CPLBDR
        END
        SUBROUTINE CPCLDM
        END
        SUBROUTINE CPLBAM
        END
        SUBROUTINE CPGETR
        END
        SUBROUTINE CPPKCL
        END
        SUBROUTINE CPGETI
        END
        SUBROUTINE CPRECT
        END
        SUBROUTINE ISHIFT
        END
        SUBROUTINE MAPTRN
        END
        SUBROUTINE PCSETR
        END
        SUBROUTINE PCGETR
        END
        SUBROUTINE PLCHHQ
        END
        SUBROUTINE PCSETI
        END
        SUBROUTINE PCGETI
        END
        SUBROUTINE SFWRLD
        END
        SUBROUTINE SFSETR
        END
        SUBROUTINE SFSETI
        END
        SUBROUTINE CURVED
        END
        SUBROUTINE GETSET
        END
        SUBROUTINE AGSETR
        END
        BLOCK DATA VELDAT
        END

ext: f


subroutine apply_linearized_bcs(r, p)

! Apply linearized boundary conditions to the Newton perturbation.

    implicit none

    integer :: mx, my, mbc, meqn
    double precision :: x_lower, y_lower, dx, dy
    common /claw_config/ mx, my, mbc, x_lower, y_lower, dx, dy, meqn

    character(len=2), dimension(4, 10) :: bc_options
    common /bc_config/ bc_options

    double precision, dimension(1-mbc:mx+mbc, 1-mbc:my+mbc, meqn), intent(in) :: r
    double precision, dimension(1-mbc:mx+mbc, 1-mbc:my+mbc, meqn), intent(inout) :: p
    integer :: i

    do i = 1, meqn
        call fill_component_ghost_cells_homogeneous(bc_options(:, i), p(:,:,i))
    end do

end subroutine apply_linearized_bcs

subroutine apply_bcs(t, q)

! Fill ghost cells with appropriate boundary values.

    implicit none

    integer :: mx, my, mbc, meqn
    double precision :: x_lower, y_lower, dx, dy
    common /claw_config/ mx, my, mbc, x_lower, y_lower, dx, dy, meqn

    character(len=2), dimension(4, 10) :: bc_options
    common /bc_config/ bc_options

    double precision, intent(in) :: t
    double precision, dimension(1-mbc:mx+mbc, 1-mbc:my+mbc, meqn), intent(inout) :: q

    double precision, dimension(2, my, meqn) :: x_lower_values, x_upper_values
    double precision, dimension(2, mx, meqn) :: y_lower_values, y_upper_values
    integer :: i

    call set_implicit_boundary_data(t, x_lower_values, x_upper_values,  &
                                    y_lower_values, y_upper_values)

    do i = 1, meqn
        call fill_component_ghost_cells(bc_options(:, i), x_lower_values(:, :, i), x_upper_values(:, :, i),  &
                              y_lower_values(:, :, i), y_upper_values(:, :, i), q(:, :, i))
    end do

end subroutine apply_bcs

subroutine fill_component_ghost_cells_homogeneous(bc_options, q_component)

! Fill ghost cells with homogeneous boundary conditions.

    implicit none

    integer :: mx, my, mbc, meqn
    double precision :: x_lower, y_lower, dx, dy
    common /claw_config/ mx, my, mbc, x_lower, y_lower, dx, dy, meqn

    character(len=2), dimension(4), intent(in) :: bc_options
    double precision, dimension(1-mbc:mx+mbc, 1-mbc:my+mbc), intent(inout) :: q_component
    double precision, dimension(:, :), allocatable :: zeros
    integer :: i

    allocate(zeros(2, max(mx, my)))
    zeros = 0.d0

    ! Dimensions are correct; there are my horizontal strips and mx vertical strips.
    call fill_component_ghost_cells(bc_options, zeros(:, 1:my), zeros(:, 1:my),  &
                          zeros(:, 1:mx), zeros(:, 1:mx), q_component)

end subroutine fill_component_ghost_cells_homogeneous

! TODO: Move to separate file, and add dedicated homogeneous method for
! convenience.
subroutine fill_component_ghost_cells(bc_options, x_lower_values, x_upper_values,  &
                        y_lower_values, y_upper_values, q)

! Fills ghost cells on either side of the domain for a single solution
! component.
!
! 'option' arguments indicate the condition used to fill the cell.
! Allowable values are:
!   'p': periodic
!   '0' or '1': Match derivative 0 or 1
!   '00', '01', '02', '03', '12', '13', or '23': Match derivatives of the
!     two specified orders.
!
! If either option is 'p', then both must be.
!
! 'value' arguments are double precision, indicating the value to be
! matched.  If the corresponding option is 'p', this value is unused.

    implicit none

    integer :: mx, my, mbc, meqn
    double precision :: x_lower, y_lower, dx, dy
    common /claw_config/ mx, my, mbc, x_lower, y_lower, dx, dy, meqn

    character(len=2), dimension(4), intent(in) :: bc_options
    double precision, intent(in), dimension(2, my) :: x_lower_values, x_upper_values
    double precision, intent(in), dimension(2, mx) :: y_lower_values, y_upper_values
    double precision, intent(inout) :: q(1-mbc:mx+mbc, 1-mbc:my+mbc)

    double precision, dimension(2, mbc) :: zeros
    double precision, dimension(2, 1) :: extrap_boundary_values
    double precision, dimension(2, 2) :: temp_corner

    zeros = 0.d0

    ! TODO: Make sure there's a safety check on the periodic option.
    call set_x_lower_ghost_cells(bc_options(1), x_lower_values, 1, my)
    call set_x_upper_ghost_cells(bc_options(2), x_upper_values, 1, my)
    call set_y_lower_ghost_cells(bc_options(3), y_lower_values, 1, mx)
    call set_y_upper_ghost_cells(bc_options(4), y_upper_values, 1, mx)

    ! Handle corners if one or both dimensions are periodic.
    if (bc_options(1) == 'p') then
        call set_x_lower_ghost_cells(bc_options(1), zeros, 1-mbc, 0)
        call set_x_lower_ghost_cells(bc_options(1), zeros, my+1, my+mbc)
        call set_x_upper_ghost_cells(bc_options(2), zeros, 1-mbc, 0)
        call set_x_upper_ghost_cells(bc_options(2), zeros, my+1, my+mbc)
    else if (bc_options(3) == 'p') then
        call set_y_lower_ghost_cells(bc_options(3), zeros, 1-mbc, 0)
        call set_y_lower_ghost_cells(bc_options(3), zeros, mx+1, mx+mbc)
        call set_y_upper_ghost_cells(bc_options(4), zeros, 1-mbc, 0)
        call set_y_upper_ghost_cells(bc_options(4), zeros, mx+1, mx+mbc)
    else
        ! x lower, y lower
        call quintic_extrap(q(1:6, -1), q(-1:0, -1))
        call quintic_extrap(q(1:6, 0), q(-1:0, 0))
        temp_corner = q(-1:0, -1:0)
        call quintic_extrap(q(-1, 1:6), q(-1, -1:0))
        call quintic_extrap(q(0, 1:6), q(0, -1:0))
        q(-1:0, -1:0) = (q(-1:0, -1:0) + temp_corner) / 2.d0

        ! x upper, y lower
        call quintic_extrap(q(mx:mx-5:-1, -1), q(mx+2:mx+1:-1, -1))
        call quintic_extrap(q(mx:mx-5:-1, 0), q(mx+2:mx+1:-1, 0))
        temp_corner = q(mx+1:mx+2, -1:0)
        call quintic_extrap(q(mx+1, 1:6), q(mx+1, -1:0))
        call quintic_extrap(q(mx+2, 1:6), q(mx+2, -1:0))
        q(mx+1:mx+2, -1:0) = (q(mx+1:mx+2, -1:0) + temp_corner) / 2.d0

        ! x lower, y upper
        call quintic_extrap(q(1:6, my+1), q(-1:0, my+1))
        call quintic_extrap(q(1:6, my+2), q(-1:0, my+2))
        temp_corner = q(-1:0, my+1:my+2)
        call quintic_extrap(q(-1, my:my-5:-1), q(-1, my+2:my+1:-1))
        call quintic_extrap(q(0, my:my-5:-1), q(0, my+2:my+1:-1))
        q(-1:0, my+1:my+2) = (q(-1:0, my+1:my+2) + temp_corner) / 2.d0

        ! x upper, y upper
        call quintic_extrap(q(mx:mx-5:-1, my+1), q(mx+2:mx+1:-1, my+1))
        call quintic_extrap(q(mx:mx-5:-1, my+2), q(mx+2:mx+1:-1, my+2))
        temp_corner = q(mx+1:mx+2, my+1:my+2)
        call quintic_extrap(q(mx+1, my:my-5:-1), q(mx+1, my+2:my+1:-1))
        call quintic_extrap(q(mx+2, my:my-5:-1), q(mx+2, my+2:my+1:-1))
        q(mx+1:mx+2, my+1:my+2) = (q(mx+1:mx+2, my+1:my+2) + temp_corner) / 2.d0
    end if

    contains

    subroutine set_x_lower_ghost_cells(options, boundary_values, iy_low, iy_high)
        implicit none
        character(len=2), intent(in) :: options
        integer, intent(in) :: iy_low, iy_high
        double precision, dimension(2, iy_low:iy_high), intent(in) :: boundary_values
        double precision, dimension(2, 2) :: boundary_coeffs
        double precision, dimension(3, 2) :: cell_coeffs
        integer :: iy

        if (options == 'p') then
            do iy = iy_low, iy_high
                q(1-mbc:0, iy) = q(mx+1-mbc:mx, iy)
            end do
        else
            call get_lower_bc_coefficients(options, dx, boundary_coeffs, cell_coeffs)
            if (len_trim(options) == 1) then
                do iy = iy_low, iy_high
                    q(-1:0, iy) = boundary_values(1, iy) * boundary_coeffs(1, :)  &
                        + matmul(q(1:3, iy), cell_coeffs)
                end do
            else
                do iy = iy_low, iy_high
                    q(-1:0, iy) = matmul(boundary_values(:, iy), boundary_coeffs)  &
                        + matmul(q(1:3, iy), cell_coeffs)
                end do
            end if
        end if
    end subroutine set_x_lower_ghost_cells

    subroutine set_x_upper_ghost_cells(options, boundary_values, iy_low, iy_high)
        implicit none
        character(len=2), intent(in) :: options
        integer, intent(in) :: iy_low, iy_high
        double precision, dimension(2, iy_low:iy_high), intent(in) :: boundary_values
        double precision, dimension(2, 2) :: boundary_coeffs
        double precision, dimension(3, 2) :: cell_coeffs
        integer :: iy

        if (options == 'p') then
            do iy = iy_low, iy_high
                q(mx+1:mx+mbc, iy) = q(1:mbc, iy)
            end do
        else
            call get_upper_bc_coefficients(options, dx, boundary_coeffs, cell_coeffs)
            if (len_trim(options) == 1) then
                do iy = iy_low, iy_high
                    q(mx+1:mx+2, iy) = boundary_values(1, iy) * boundary_coeffs(1, :)  &
                        + matmul(q(mx-2:mx, iy), cell_coeffs)
                end do
            else
                do iy = iy_low, iy_high
                    q(mx+1:mx+2, iy) = matmul(boundary_values(:, iy), boundary_coeffs)  &
                        + matmul(q(mx-2:mx, iy), cell_coeffs)
                end do
            end if
        end if
    end subroutine set_x_upper_ghost_cells

    subroutine set_y_lower_ghost_cells(options, boundary_values, ix_low, ix_high)
        implicit none
        character(len=2), intent(in) :: options
        integer, intent(in) :: ix_low, ix_high
        double precision, dimension(2, ix_low:ix_high), intent(in) :: boundary_values
        double precision, dimension(2, 2) :: boundary_coeffs
        double precision, dimension(3, 2) :: cell_coeffs
        integer :: ix

        if (options == 'p') then
            do ix = ix_low, ix_high
                q(ix, 1-mbc:0) = q(ix, my+1-mbc:my)
            end do
        else
            call get_lower_bc_coefficients(bc_options(3), dy, boundary_coeffs, cell_coeffs)
            if (len_trim(options) == 1) then
                do ix = ix_low, ix_high
                    q(ix, -1:0) = boundary_values(1, ix) * boundary_coeffs(1, :)  &
                        + matmul(q(ix, 1:3), cell_coeffs)
                end do
            else
                do ix = ix_low, ix_high
                    q(ix, -1:0) = matmul(boundary_values(:, ix), boundary_coeffs)  &
                        + matmul(q(ix, 1:3), cell_coeffs)
                end do
            end if
        end if
    end subroutine set_y_lower_ghost_cells

    subroutine set_y_upper_ghost_cells(options, boundary_values, ix_low, ix_high)
        implicit none
        character(len=2), intent(in) :: options
        integer, intent(in) :: ix_low, ix_high
        double precision, dimension(2, ix_low:ix_high), intent(in) :: boundary_values
        double precision, dimension(2, 2) :: boundary_coeffs
        double precision, dimension(3, 2) :: cell_coeffs
        integer :: ix

        if (options == 'p') then
            do ix = ix_low, ix_high
                q(ix, my+1:my+mbc) = q(ix, 1:mbc)
            end do
        else
            call get_upper_bc_coefficients(bc_options(4), dy, boundary_coeffs, cell_coeffs)
            if (len_trim(options) == 1) then
                do ix = ix_low, ix_high
                    q(ix, my+1:my+2) = boundary_values(1, ix) * boundary_coeffs(1, :)  &
                        + matmul(q(ix, my-2:my), cell_coeffs)
                end do
            else
                do ix = 1, mx
                    q(ix, my+1:my+2) = matmul(boundary_values(:, ix), boundary_coeffs)  &
                        + matmul(q(ix, my-2:my), cell_coeffs)
                end do
            end if
        end if
    end subroutine set_y_upper_ghost_cells

end subroutine fill_component_ghost_cells

subroutine quintic_extrap(source, dest)
    implicit none
    double precision, dimension(6), intent(in) :: source
    double precision, dimension(2), intent(out) :: dest
    dest(1) = 21.d0 * source(1) - 70.d0 * source(2) + 105.d0 * source(3)  &
        - 84.d0 * source(4) + 35.d0 * source(5) - 6.d0 * source(6)
    dest(2) = 6.d0 * source(1) - 15.d0 * source(2) + 20.d0 * source(3)  &
        - 15.d0 * source(4) + 6.d0 * source(5) - source(6)
end subroutine quintic_extrap


text:
!Copyright (c) 2009-2016 Tom Schoonjans
!All rights reserved.

!Redistribution and use in source and binary forms, with or without
!modification, are permitted provided that the following conditions are met:
!    * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
!    * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
!    * The names of the contributors may not be used to endorse or promote products derived from this software without specific prior written permission.

!THIS SOFTWARE IS PROVIDED BY Tom Schoonjans ''AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL Tom Schoonjans BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

PROGRAM xrltest

USE xraylib

IMPLICIT NONE

TYPE (compoundData), POINTER :: cd

CHARACTER (KIND=C_CHAR,LEN=10) :: compound1 = 'Ca(HCO3)2'
CHARACTER (KIND=C_CHAR,LEN=5) :: compound2 = 'SiO2'
INTEGER :: i
TYPE (Crystal_Struct), POINTER :: cryst
REAL (C_DOUBLE) :: bragg, q, energy, debye_temp_factor, f0, fp, fpp,&
rel_angle
REAL (C_DOUBLE) :: dw
COMPLEX (C_DOUBLE) :: F_H, F_0, F_Hbar
REAL (C_DOUBLE), PARAMETER :: PI = 4.D0*DATAN(1.D0)

TYPE (compoundDataNIST), POINTER :: cdn
TYPE (radioNuclideData), POINTER :: rnd
CHARACTER (KIND=C_CHAR, LEN=NIST_LIST_STRING_LENGTH), POINTER, &
DIMENSION(:) :: nistCompounds
CHARACTER (KIND=C_CHAR, LEN=RADIO_NUCLIDE_STRING_LENGTH), POINTER, &
DIMENSION(:) :: radioNuclides
CHARACTER (KIND=C_CHAR, LEN=CRYSTAL_STRING_LENGTH), POINTER, &
DIMENSION(:) :: crystals

CALL XRayInit()

WRITE (6,'(A)') 'Example of fortran program using xraylib'
WRITE (6,'(A,F12.6,A)') 'Density of pure Al: ',ElementDensity(13),' g/cm3'
WRITE (6,'(A,F12.6)') 'Ca K-alpha Fluorescence Line Energy: ',LineEnergy(20,KA_LINE)
WRITE (6,'(A,F12.6)') 'Fe partial photoionization cs of L3 at 6.0 keV: ',CS_Photo_Partial(26,L3_SHELL,6.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Zr L1 edge energy: ',EdgeEnergy(40,L1_SHELL)
WRITE (6,'(A,F12.6)') 'Pb Lalpha XRF production cs at 20.0 keV (jump approx): ',CS_FluorLine(82,LA_LINE,20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb Lalpha XRF production cs at 20.0 keV (Kissel): ',CS_FluorLine_Kissel(82,LA_LINE,20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Bi M1N2 radiative rate: ',RadRate(83,M1N2_LINE)
WRITE (6,'(A,F12.6)') 'U M3O3 Fluorescence Line Energy: ',LineEnergy(92,M3O3_LINE)

!CompoundParser tests
cd => CompoundParser(compound1)
IF (.NOT. ASSOCIATED(cd)) THEN
        CALL EXIT(1)
ENDIF
WRITE (6,'(A,F12.6,A,I4,A,F12.6,A)') 'Ca(HCO3)2 contains ',cd%nAtomsAll,' atoms, ',&
cd%nElements,' elements and has a molar mass of ', cd%molarMass, ' g/mol'
DO i=1,cd%nElements
        WRITE (6,'(A,I2,A,F12.6,A,F12.6,A)') 'Element ',cd%Elements(i),' : ',&
        cd%massFractions(i)*100.0_C_DOUBLE,' % and ',&
        cd%nAtoms(i)
ENDDO

!Free the memory allocated for the arrays
DEALLOCATE(cd)

cd => CompoundParser(compound2)
IF (.NOT. ASSOCIATED(cd)) THEN
        CALL EXIT(1)
ENDIF
WRITE (6,'(A,F12.6,A,I4,A,F12.6,A)') 'SiO2 contains ',cd%nAtomsAll,' atoms, ',&
cd%nElements,' elements and has a molar mass of ', cd%molarMass, ' g/mol'
DO i=1,cd%nElements
        WRITE (6,'(A,I2,A,F12.6,A,F12.6,A)') 'Element ',cd%Elements(i),' : ',&
        cd%massFractions(i)*100.0_C_DOUBLE,' % and ',&
        cd%nAtoms(i)
ENDDO

!Free the memory allocated for the arrays
DEALLOCATE(cd)

WRITE (6,'(A,F12.6)') 'Ca(HCO3)2 Rayleigh cs at 10.0 keV: ',CS_Rayl_CP('Ca(HCO3)2',10.0_C_DOUBLE)

WRITE (6,'(A,ES14.6,A,ES14.6,A)') 'CS2 Refractive Index at 10.0 keV : ', &
        Refractive_Index_Re('CS2',10.0_C_DOUBLE,1.261_C_DOUBLE),' - ',&
        Refractive_Index_Im('CS2',10.0_C_DOUBLE,1.261_C_DOUBLE),' i'
WRITE (6,'(A,ES14.6,A,ES14.6,A)') 'C16H14O3 Refractive Index at 1 keV : ', &
        Refractive_Index_Re('C16H14O3',1.0_C_DOUBLE,1.2_C_DOUBLE),' - ',&
        Refractive_Index_Im('C16H14O3',1.0_C_DOUBLE,1.2_C_DOUBLE),' i'
WRITE (6,'(A,ES14.6,A,ES14.6,A)') 'SiO2 Refractive Index at 5.0 keV : ', &
        Refractive_Index_Re('SiO2',5.0_C_DOUBLE,2.65_C_DOUBLE),' - ',&
        Refractive_Index_Im('SiO2',5.0_C_DOUBLE,2.65_C_DOUBLE),' i'
WRITE (6,'(A,F12.6)') 'Compton profile for Fe at pz = 1.1 : ' ,&
        ComptonProfile(26,1.1_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'M5 Compton profile for Fe at pz = 1.1 : ' ,&
        ComptonProfile_Partial(26,M5_SHELL,1.1_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'K atomic level width for Fe: ',&
        AtomicLevelWidth(26,K_SHELL)
WRITE (6,'(A,F12.6)') 'Bi L2-M5M5 Auger non-radiative rate: ',&
        AugerRate(86,L2_M5M5_AUGER)
WRITE (6,'(A,F12.6)') 'Bi L3 Auger yield: ',&
        AugerYield(86, L3_SHELL)

WRITE (6,'(A,F12.6)') 'M1->M5 Coster-Kronig transition probability for Au : ',CosKronTransProb(79,FM15_TRANS)
WRITE (6,'(A,F12.6)') 'L1->L3 Coster-Kronig transition probability for Fe : ',CosKronTransProb(26,FL13_TRANS)
WRITE (6,'(A,F12.6)') 'Au Ma1 XRF production cs at 10.0 keV (Kissel): ',CS_FluorLine_Kissel(79,MA1_LINE,10.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Au Mb XRF production cs at 10.0 keV (Kissel): ',CS_FluorLine_Kissel(79,MB_LINE,10.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Au Mg XRF production cs at 10.0 keV (Kissel): ',CS_FluorLine_Kissel(79,MG_LINE,10.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb Malpha XRF production cs at 20.0 keV with cascade effect: ',&
CS_FluorLine_Kissel(82,MA1_LINE,20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb Malpha XRF production cs at 20.0 keV with radiative cascade effect: ',&
CS_FluorLine_Kissel_Radiative_Cascade(82,MA1_LINE,20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb Malpha XRF production cs at 20.0 keV with non-radiative cascade effect: ',&
CS_FluorLine_Kissel_Nonradiative_Cascade(82,MA1_LINE,20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb Malpha XRF production cs at 20.0 keV without cascade effect: ',&
CS_FluorLine_Kissel_no_Cascade(82,MA1_LINE,20.0_C_DOUBLE)

WRITE (6,'(A,F12.6)') 'Al mass energy-absorption cs at 20.0 keV: ', CS_Energy(13, 20.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'Pb mass energy-absorption cs at 40.0 keV: ', CS_Energy(82, 40.0_C_DOUBLE)
WRITE (6,'(A,F12.6)') 'CdTe mass energy-absorption cs at 40.0 keV: ',&
CS_Energy_CP('CdTe', 40.0_C_DOUBLE)

WRITE (6, '(A,F12.6)') 'Sr anomalous scattering factor Fi at 10.0 keV: ', Fi(38, 10.0_C_DOUBLE)
WRITE (6, '(A,F12.6)') 'Sr anomalous scattering factor Fii at 10.0 keV: ', Fii(38, 10.0_C_DOUBLE)

WRITE (6,'(A,A)') 'Symbol of element 26 is: ',AtomicNumberToSymbol(26)
WRITE (6,'(A,I3)') 'Number of element Fe is: ',SymbolToAtomicNumber('Fe')

cryst => Crystal_GetCrystal('Si')
IF (.NOT.ASSOCIATED(cryst)) CALL EXIT(1)

WRITE (6,'(A,3F12.3)') 'Si unit cell dimensions are ',cryst%a,cryst%b,cryst%c
WRITE (6,'(A,3F12.3)') 'Si unit cell angles are ',&
cryst%alpha,cryst%beta,cryst%gamma
WRITE (6,'(A,F12.3)') 'Si unit cell volume is ',cryst%volume
WRITE (6,'(A)') 'Si atoms at '
WRITE (6,'(A)') ' Z  fraction    X        Y        Z'
DO i=1,cryst%n_atom
        WRITE (6, '(I3,4F9.3)') cryst%atom(i)%Zatom, cryst%atom(i)%fraction,&
        cryst%atom(i)%x, cryst%atom(i)%y, cryst%atom(i)%z
ENDDO

WRITE (6,'(A)') ''

!Si diffraction parameters
WRITE (6, '(A)') 'Si111 at 8 KeV. Incidence at the Bragg angle:'
bragg = Bragg_angle(cryst, 8.0_C_DOUBLE, 1, 1, 1)
WRITE (6, '(A,F12.6,A,F12.6)') '  Bragg angle: Rad: ',bragg,' Deg: ',&
bragg*180.0/PI
q = Q_scattering_amplitude (cryst, 8.0_C_DOUBLE, 1, 1, 1, 1.0_C_DOUBLE)
WRITE (6, '(A, F12.6)') '  Q Scattering amplitude: ',q
energy = 8.0
debye_temp_factor = 1.0
i = Atomic_Factors (14, energy, q, debye_temp_factor, f0, fp, fpp)

WRITE (6, '(A,F12.6,A,F12.6,A,F12.6)')&
'  Atomic factors (Z=14) f0, fp, fpp: ', f0, ', ',fp, ', i*',fpp

rel_angle = 1.0

F_H = Crystal_F_H_StructureFactor (cryst, energy, 1, 1, 1, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  FH(1,1,1) structure factor: (',&
REAL(F_H),', ',AIMAG(F_H),')'

F_0 = Crystal_F_H_StructureFactor (cryst, energy, 0, 0, 0, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  F0=FH(0,0,0) structure factor: (',&
REAL(F_0),', ',AIMAG(F_0),')'

DEALLOCATE(cryst)

WRITE (6,'(A)') ''

! Diamond diffraction parameters
cryst => Crystal_GetCrystal('Diamond')
IF (.NOT.ASSOCIATED(cryst)) CALL EXIT(1)
WRITE (6,'(A)') 'Diamond 111 at 8 KeV. Incidence at the Bragg angle:'

bragg = Bragg_angle (cryst, energy, 1, 1, 1)
WRITE (6, '(A,F12.6,A,F12.6)') '  Bragg angle: Rad: ',bragg,' Deg: ',&
bragg*180.0/PI

q = Q_scattering_amplitude (cryst, 8.0_C_DOUBLE, 1, 1, 1, rel_angle)
WRITE (6, '(A, F12.6)') '  Q Scattering amplitude: ',q
energy = 8.0
debye_temp_factor = 1.0
i = Atomic_Factors (6, energy, q, debye_temp_factor, f0, fp, fpp)

WRITE (6, '(A,F12.6,A,F12.6,A,F12.6)')&
'  Atomic factors (Z=6) f0, fp, fpp: ', f0, ', ',fp, ', i*',fpp

F_H = Crystal_F_H_StructureFactor (cryst, energy, 1, 1, 1, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  FH(1,1,1) structure factor: (',&
REAL(F_H),', ',AIMAG(F_H),')'

F_0 = Crystal_F_H_StructureFactor (cryst, energy, 0, 0, 0, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  F0=FH(0,0,0) structure factor: (',&
REAL(F_0),', ',AIMAG(F_0),')'

F_Hbar = Crystal_F_H_StructureFactor (cryst, energy, -1, -1, -1, debye_temp_factor,&
rel_angle)

dw = 1e10 * 2 * (R_E / cryst%volume) * &
(KEV2ANGST * KEV2ANGST/ (energy *energy)) * &
SQRT(ABS(F_H * F_Hbar)) / PI / SIN(2*bragg)

WRITE (6, '(A,F12.6,A)') '  Darwin width: ',1.0E6*dw,' micro-radians'

DEALLOCATE(cryst)

WRITE (6,'(A)') ''

! Alpha Quartz diffraction parameters

cryst => Crystal_GetCrystal('AlphaQuartz')
IF (.NOT.ASSOCIATED(cryst)) CALL EXIT(1)
WRITE (6, '(A)') 'Alpha Quartz 020 at 8 KeV. Incidence at the Bragg angle:'

bragg = Bragg_angle (cryst, energy, 0, 2, 0)
WRITE (6, '(A,F12.6,A,F12.6)') '  Bragg angle: Rad: ',bragg,' Deg: ',&
bragg*180.0/PI

q = Q_scattering_amplitude (cryst, 8.0_C_DOUBLE, 0, 2, 0, rel_angle)
WRITE (6, '(A, F12.6)') '  Q Scattering amplitude: ',q
i = Atomic_Factors (8, energy, q, debye_temp_factor, f0, fp, fpp)

WRITE (6, '(A,F12.6,A,F12.6,A,F12.6)')&
'  Atomic factors (Z=8) f0, fp, fpp: ', f0, ', ',fp, ', i*',fpp

F_H = Crystal_F_H_StructureFactor (cryst, energy, 0, 2, 0, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  FH(0,2,0) structure factor: (',&
REAL(F_H),', ',AIMAG(F_H),')'

F_0 = Crystal_F_H_StructureFactor (cryst, energy, 0, 0, 0, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  F0=FH(0,0,0) structure factor: (',&
REAL(F_0),', ',AIMAG(F_0),')'

WRITE (6,'(A)') ''

! Muscovite diffraction parameters

cryst => Crystal_GetCrystal('Muscovite')
IF (.NOT.ASSOCIATED(cryst)) CALL EXIT(1)
WRITE (6, '(A)') 'Muscovite 331 at 8 KeV. Incidence at the Bragg angle:'

bragg = Bragg_angle (cryst, energy, 3, 3, 1)
WRITE (6, '(A,F12.6,A,F12.6)') '  Bragg angle: Rad: ',bragg,' Deg: ',&
bragg*180.0/PI

q = Q_scattering_amplitude (cryst, 8.0_C_DOUBLE, 3, 3, 1, rel_angle)
WRITE (6, '(A, F12.6)') '  Q Scattering amplitude: ',q
i = Atomic_Factors (19, energy, q, debye_temp_factor, f0, fp, fpp)

WRITE (6, '(A,F12.6,A,F12.6,A,F12.6)')&
'  Atomic factors (Z=19) f0, fp, fpp: ', f0, ', ',fp, ', i*',fpp

F_H = Crystal_F_H_StructureFactor (cryst, energy, 3, 3, 1, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  FH(3,3,1) structure factor: (',&
REAL(F_H),', ',AIMAG(F_H),')'

F_0 = Crystal_F_H_StructureFactor (cryst, energy, 0, 0, 0, debye_temp_factor,&
rel_angle)
WRITE (6, '(A,F12.6,A,F12.6,A)') '  F0=FH(0,0,0) structure factor: (',&
REAL(F_0),', ',AIMAG(F_0),')'

DEALLOCATE(cryst)

crystals => Crystal_GetCrystalsList()
WRITE (6, '(A)') 'List of available crystals'
DO i=1,SIZE(crystals)
        WRITE (6, '(A,I3,A,A)') '  Crystal',i,': ',TRIM(crystals(i))
ENDDO

WRITE (6, '(A)') ''

cdn => GetCompoundDataNISTByName('Uranium Monocarbide')
WRITE (6, '(A)') 'Uranium Monocarbide'
WRITE (6, '(A,A)') '  Name: ', TRIM(cdn%name)
WRITE (6, '(A,F12.6,A)') '  Density: ',cdn%density,' g/cm3'
DO i=1,cdn%nElements
        WRITE (6, '(A,I2,A,F12.6,A)') '  Element ',&
        cdn%Elements(i),': ', &
        cdn%massFractions(i)*100.0, ' %'
ENDDO
DEALLOCATE(cdn)

cdn => GetCompoundDataNISTByIndex(NIST_COMPOUND_BRAIN_ICRP)
WRITE (6, '(A)') 'NIST_COMPOUND_BRAIN_ICRP'
WRITE (6, '(A,A)') '  Name: ', TRIM(cdn%name)
WRITE (6, '(A,F12.6,A)') '  Density: ',cdn%density,' g/cm3'
DO i=1,cdn%nElements
        WRITE (6, '(A,I2,A,F12.6,A)') '  Element ',&
        cdn%Elements(i),': ', &
        cdn%massFractions(i)*100.0, ' %'
ENDDO
DEALLOCATE(cdn)

nistCompounds => GetCompoundDataNISTList()
WRITE (6, '(A)') 'List of available NIST compounds'
DO i=1,SIZE(nistCompounds)
        WRITE (6, '(A,I3,A,A)') '  Compound ',i,': ',TRIM(nistCompounds(i))
ENDDO

DEALLOCATE(nistCompounds)

WRITE (6, '(A)') ''

rnd => GetRadioNuclideDataByName('109Cd')
WRITE (6, '(A)') '109Cd'
WRITE (6, '(A,A)') '  Name: ', TRIM(rnd%name)
WRITE (6, '(A,I3)') '  Z: ',rnd%Z
WRITE (6, '(A,I3)') '  A: ',rnd%A
WRITE (6, '(A,I3)') '  N: ',rnd%N
WRITE (6, '(A,I3)') '  Z_xray: ',rnd%Z_xray
WRITE (6, '(A)') '  X-rays: '
DO i=1,rnd%nXrays
        WRITE (6, '(A,F12.6,A,F12.6)') '  ',&
        LineEnergy(rnd%Z_xray, rnd%XrayLines(i)),' keV -> ', &
        rnd%XrayIntensities(i)
ENDDO
WRITE (6, '(A)') '  Gamma rays: '
DO i=1,rnd%nGammas
        WRITE (6, '(A,F12.6,A,F12.6)') '  ',&
        rnd%GammaEnergies(i),' keV -> ', &
        rnd%GammaIntensities(i)
ENDDO
DEALLOCATE(rnd)

rnd => GetRadioNuclideDataByIndex(RADIO_NUCLIDE_125I)
WRITE (6, '(A)') 'RADIO_NUCLIDE_125I'
WRITE (6, '(A,A)') '  Name: ', TRIM(rnd%name)
WRITE (6, '(A,I3)') '  Z: ',rnd%Z
WRITE (6, '(A,I3)') '  A: ',rnd%A
WRITE (6, '(A,I3)') '  N: ',rnd%N
WRITE (6, '(A,I3)') '  Z_xray: ',rnd%Z_xray
WRITE (6, '(A)') '  X-rays: '
DO i=1,rnd%nXrays
        WRITE (6, '(A,F12.6,A,F12.6)') '  ',&
        LineEnergy(rnd%Z_xray, rnd%XrayLines(i)),' keV -> ', &
        rnd%XrayIntensities(i)
ENDDO
WRITE (6, '(A)') '  Gamma rays: '
DO i=1,rnd%nGammas
        WRITE (6, '(A,F12.6,A,F12.6)') '  ',&
        rnd%GammaEnergies(i),' keV -> ', &
        rnd%GammaIntensities(i)
ENDDO
DEALLOCATE(rnd)

radioNuclides => GetRadioNuclideDataList()
WRITE (6, '(A)') 'List of available radionuclides'
DO i=1,SIZE(radioNuclides)
        WRITE (6, '(A,I3,A,A)') '  Radionuclide',i,': ',TRIM(radioNuclides(i))
ENDDO

DEALLOCATE(radioNuclides)
WRITE (6,'(A)') ''
WRITE (6,'(A)') '--------------------------- END OF XRLEXAMPLE3 -------------------------------'
WRITE (6,'(A)') ''
ENDPROGRAM

max_stars_repo_name: mschollmeier/xraylib

max_issues_repo_name: mschollmeier/xraylib

max_forks_repo_name: mschollmeier/xraylib


text: #redirect wiki:sacramento:OsakaYa

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      program bucle_simultaneo
        implicit none
        real, parameter :: pi =3.14159265
        integer, parameter :: n=10
        real :: resultado_sin(n)
        integer :: i
        do  concurrent (i=1:n) !Cuidado la sintaxis es ligeramente diferente
          resultado_sin(i)=sin(i*pi/4.)
        end do
        print *, resultado_sin
      end program bucle_simultaneo

max_stars_repo_name: Marc-xyz/InicioRapidoEnFortran

max_issues_repo_name: Marc-xyz/InicioRapidoEnFortran

max_forks_repo_name: Marc-xyz/InicioRapidoEnFortran


text:
      SUBROUTINE ALG18 (BETA1,BETA2,I1,I2,FACT,X0,Y0,S0,XR,Y1,X1,Y2,RDI
     1US,S,C1)
C
      DIMENSION S(80)
C
      DELX=XR/FLOAT(I2-I1)
      XX=X0
      I3=I1+1
      IF (BETA1.EQ.BETA2) GO TO 20
      Y1=-COS(BETA1/C1)/(SIN(BETA1/C1)-SIN(BETA2/C1))
      X1=SIN(BETA1/C1)/(SIN(BETA1/C1)-SIN(BETA2/C1))
      Y2=TAN((BETA1+BETA2)/(2.0*C1))
      RDIUS=ABS(1.0/(SIN(BETA1/C1)-SIN(BETA2/C1)))
      Y2=Y2*FACT+Y0
      Y1=Y1*FACT+Y0
      X1=X1*FACT+X0
      RDIUS=RDIUS*FACT
      DO 10 J=I3,I2
      XX=XX+DELX
      PHI1=ATAN(-1./SQRT(RDIUS**2-(XX-X1)**2)*(XX-X1))
      IF ((BETA1-BETA2).LT.0.0) PHI1=-PHI1
      PHI2=ABS(BETA1/C1-PHI1)
10    S(J)=RDIUS*PHI2+S0
      RETURN
20    AM=TAN(BETA1/C1)
      DO 30 J=I3,I2
      XX=XX+DELX
30    S(J)=(XX-X0)*SQRT(AM*AM+1.0)+S0
      Y2=AM*(XX-X0)+Y0
      RETURN
      END