

























































































































 


































































































































































 






















































































































































































 






















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                





































module adm_metric
  use constants
  use tensor
  implicit none
  private
  public calc_4metric
  public calc_4metricderivs
  public calc_4metricderivs2
  public calc_3metric
  public calc_3metricderivs
  
  public calc_extcurv
  public calc_3metricdot
  public calc_3metricderivdot


  

  
contains
  
  subroutine calc_4metric (gg, alfa, beta, g4)
    REAL*8, intent(in)  :: gg(3,3), alfa, beta(3)
    REAL*8, intent(out) :: g4(0:3,0:3)
    REAL*8 :: betal(3)
    
    
    betal = matmul(gg, beta)
    g4(0  ,0  ) = -alfa**2 + sum(betal*beta)
    g4(1:3,0  ) = betal
    g4(0  ,1:3) = betal
    g4(1:3,1:3) = gg
  end subroutine calc_4metric
  
  subroutine calc_4metricderivs (gg,alfa,beta, dgg,dalfa,dbeta, &
       gg_dot,alfa_dot,beta_dot, g4,dg4)
    REAL*8, intent(in)  :: gg(3,3),alfa,beta(3)
    REAL*8, intent(in)  :: dgg(3,3,3),dalfa(3),dbeta(3,3)
    REAL*8, intent(in)  :: gg_dot(3,3),alfa_dot,beta_dot(3)
    REAL*8, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
    REAL*8 :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3)
    REAL*8 :: betal(3),d4betal(3,0:3)
    integer   :: i,j
    
    
    forall (i=1:3)
       betal(i) = sum(gg(i,:) * beta(:))
    end forall
    g4(0  ,0  ) = -alfa**2 + sum(betal(:) * beta(:))
    g4(1:3,0  ) = betal
    g4(0  ,1:3) = betal
    g4(1:3,1:3) = gg
    
    
    d4gg  (:,:,0  ) = gg_dot(:,:)
    d4gg  (:,:,1:3) = dgg(:,:,:)
    d4alfa(    0  ) = alfa_dot
    d4alfa(    1:3) = dalfa(:)
    d4beta(:,  0  ) = beta_dot(:)
    d4beta(:,  1:3) = dbeta(:,:)
    
    forall (i=1:3, j=0:3)
       d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j))
    end forall
    
    forall (i=0:3)
       dg4(0  ,0  ,i) = - 2 * alfa * d4alfa(i) &
            &           + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i))
       dg4(1:3,0  ,i) = d4betal(:,i)
       dg4(0  ,1:3,i) = d4betal(:,i)
       dg4(1:3,1:3,i) = d4gg(:,:,i)
    end forall
  end subroutine calc_4metricderivs
  
  subroutine calc_4metricderivs2 (gg,alfa,beta, dgg,dalfa,dbeta, &
       ddgg,ddalfa,ddbeta, gg_dot,alfa_dot,beta_dot, &
       gg_dot2,alfa_dot2,beta_dot2, dgg_dot,dalfa_dot,dbeta_dot, g4,dg4,ddg4)
    REAL*8, intent(in)  :: gg(3,3),alfa,beta(3)
    REAL*8, intent(in)  :: dgg(3,3,3),dalfa(3),dbeta(3,3)
    REAL*8, intent(in)  :: ddgg(3,3,3,3),ddalfa(3,3),ddbeta(3,3,3)
    REAL*8, intent(in)  :: gg_dot(3,3),alfa_dot,beta_dot(3)
    REAL*8, intent(in)  :: gg_dot2(3,3),alfa_dot2,beta_dot2(3)
    REAL*8, intent(in)  :: dgg_dot(3,3,3),dalfa_dot(3),dbeta_dot(3,3)
    REAL*8, intent(out) :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
    REAL*8, intent(out) :: ddg4(0:3,0:3,0:3,0:3)
    REAL*8 :: d4gg(3,3,0:3),d4alfa(0:3),d4beta(3,0:3)
    REAL*8 :: dd4gg(3,3,0:3,0:3),dd4alfa(0:3,0:3),dd4beta(3,0:3,0:3)
    REAL*8 :: betal(3),d4betal(3,0:3),dd4betal(3,0:3,0:3)
    integer   :: i,j,k
    
    
    forall (i=1:3)
       betal(i) = sum(gg(i,:) * beta(:))
    end forall
    g4(0  ,0  ) = -alfa**2 + sum(betal(:) * beta(:))
    g4(1:3,0  ) = betal
    g4(0  ,1:3) = betal
    g4(1:3,1:3) = gg
    
    
    d4gg  (:,:,0  ) = gg_dot(:,:)
    d4gg  (:,:,1:3) = dgg(:,:,:)
    d4alfa(    0  ) = alfa_dot
    d4alfa(    1:3) = dalfa(:)
    d4beta(:,  0  ) = beta_dot(:)
    d4beta(:,  1:3) = dbeta(:,:)
    
    forall (i=1:3, j=0:3)
       d4betal(i,j) = sum(d4gg(i,:,j) * beta(:) + gg(i,:) * d4beta(:,j))
    end forall
    
    forall (i=0:3)
       dg4(0  ,0  ,i) = - 2 * alfa * d4alfa(i) &
            &           + sum(d4betal(:,i) * beta(:) + betal(:) * d4beta(:,i))
       dg4(1:3,0  ,i) = d4betal(:,i)
       dg4(0  ,1:3,i) = d4betal(:,i)
       dg4(1:3,1:3,i) = d4gg(:,:,i)
    end forall
    
    
    dd4gg  (:,:,0  ,0  ) = gg_dot2(:,:)
    dd4gg  (:,:,1:3,0  ) = dgg_dot(:,:,:)
    dd4gg  (:,:,0  ,1:3) = dgg_dot(:,:,:)
    dd4gg  (:,:,1:3,1:3) = ddgg(:,:,:,:)
    dd4alfa(    0  ,0  ) = alfa_dot2
    dd4alfa(    1:3,0  ) = dalfa_dot(:)
    dd4alfa(    0  ,1:3) = dalfa_dot(:)
    dd4alfa(    1:3,1:3) = ddalfa(:,:)
    dd4beta(:,  0  ,0  ) = beta_dot2(:)
    dd4beta(:,  1:3,0  ) = dbeta_dot(:,:)
    dd4beta(:,  0  ,1:3) = dbeta_dot(:,:)
    dd4beta(:,  1:3,1:3) = ddbeta(:,:,:)
    
    
    
    
    
    forall (i=1:3, j=0:3, k=0:3)
       dd4betal(i,j,k) = sum(+ dd4gg(i,:,j,k) *    beta(:)      &
            &                +  d4gg(i,:,j)   *  d4beta(:,k)    &
            &                +  d4gg(i,:,k)   *  d4beta(:,j)    &
            &                +    gg(i,:)     * dd4beta(:,j,k))
    end forall
    
    
    
    
    
    
    
    
    
    
    
    forall (i=0:3, j=0:3)
       ddg4(0  ,0  ,i,j) = - 2 * d4alfa(j) * d4alfa(i)        &
            &              - 2 * alfa * dd4alfa(i,j)          &
            &              + sum(+ dd4betal(:,i,j) * beta(:)  &
            &                    + d4betal(:,i) * d4beta(:,j) &
            &                    + d4betal(:,j) * d4beta(:,i) &
            &                    + betal(:) * dd4beta(:,i,j))
       ddg4(1:3,0  ,i,j) = dd4betal(:,i,j)
       ddg4(0  ,1:3,i,j) = dd4betal(:,i,j)
       ddg4(1:3,1:3,i,j) = dd4gg(:,:,i,j)
    end forall
  end subroutine calc_4metricderivs2
  
  
  
  subroutine calc_3metric (g4, gg, alfa, beta)
    REAL*8, intent(in)  :: g4(0:3,0:3)
    REAL*8, intent(out) :: gg(3,3), alfa, beta(3)
    REAL*8 :: betal(3)
    REAL*8 :: dtg, gu(3,3)
    
    
    betal = g4(1:3,0)
    gg = g4(1:3,1:3)
    call calc_det (gg, dtg)
    call calc_inv (gg, dtg, gu)
    beta = matmul(gu, betal)
    alfa = sqrt(sum(betal*beta) - g4(0,0))
  end subroutine calc_3metric
  
  subroutine calc_3metricderivs (g4,dg4, gg,alfa,beta, dgg,dalfa,dbeta, &
       gg_dot,alfa_dot,beta_dot)
    REAL*8, intent(in)  :: g4(0:3,0:3),dg4(0:3,0:3,0:3)
    REAL*8, intent(out) :: gg(3,3),alfa,beta(3)
    REAL*8, intent(out) :: dgg(3,3,3),dalfa(3),dbeta(3,3)
    REAL*8, intent(out) :: gg_dot(3,3),alfa_dot,beta_dot(3)
    REAL*8 :: betal(3),d4betal(3,0:3)
    REAL*8 :: dtg,gu(3,3),dgu(3,3,3),gu_dot(3,3)
    REAL*8 :: d4gg(3,3,0:3),d4gu(3,3,0:3)
    REAL*8 :: d4alfa(0:3),d4beta(3,0:3)
    integer   :: i,j
    
    
    gg = g4(1:3,1:3)
    call calc_det (gg, dtg)
    call calc_inv (gg, dtg, gu)
    betal = g4(1:3,0)
    beta = matmul(gu, betal)
    alfa = sqrt(sum(betal*beta) - g4(0,0))
    
    forall (i=0:3)
       d4gg(:,:,i) = dg4(1:3,1:3,i)
    end forall
    gg_dot = d4gg(:,:,0)
    dgg = d4gg(:,:,1:3)
    
    call calc_invderiv (gu, dgg, dgu)
    call calc_invdot (gu, gg_dot, gu_dot)
    d4gu(:,:,0) = gu_dot
    d4gu(:,:,1:3) = dgu
    
    forall (i=0:3)
       d4betal(:,i) = dg4(1:3,0,i)
    end forall
    forall (i=0:3, j=1:3)
       d4beta(j,i) = sum(d4gu(j,:,i)*betal) + sum(gu(j,:)*d4betal(:,i))
    end forall
    forall (i=0:3)
       d4alfa(i) = 1/(2*alfa) &
            * (sum(d4betal(:,i)*beta) + sum(betal*d4beta(:,i)) - dg4(0,0,i))
    end forall
    
    alfa_dot = d4alfa(0)
    dalfa = d4alfa(1:3)
    beta_dot = d4beta(:,0)
    dbeta = d4beta(:,1:3)
  end subroutine calc_3metricderivs
  
  
  
  subroutine calc_extcurv (gg, dgg, gg_dot, alfa, beta, dbeta, kk)
    REAL*8, intent(in)  :: gg(3,3), dgg(3,3,3), gg_dot(3,3)
    REAL*8, intent(in)  :: alfa
    REAL*8, intent(in)  :: beta(3), dbeta(3,3)
    REAL*8, intent(out) :: kk(3,3)
    integer :: i,j,k
    
    
    
    do i=1,3
       do j=1,3
          kk(i,j) = - gg_dot(i,j)
          do k=1,3
             kk(i,j) = kk(i,j) &
                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
                  + beta(k) * dgg(i,j,k)
          end do
          kk(i,j) = kk(i,j) / (2*alfa)
       end do
    end do
    
  end subroutine calc_extcurv
  
  subroutine calc_3metricdot (gg, dgg, kk, alfa, beta, dbeta, gg_dot)
    REAL*8, intent(in)  :: gg(3,3), dgg(3,3,3)
    REAL*8, intent(in)  :: kk(3,3)
    REAL*8, intent(in)  :: alfa
    REAL*8, intent(in)  :: beta(3), dbeta(3,3)
    REAL*8, intent(out) :: gg_dot(3,3)
    integer :: i,j,k
    
    
    
    do i=1,3
       do j=1,3
          gg_dot(i,j) = -2*alfa * kk(i,j)
          do k=1,3
             gg_dot(i,j) = gg_dot(i,j) &
                  + gg(k,j) * dbeta(k,i) + gg(i,k) * dbeta(k,j) &
                  + beta(k) * dgg(i,j,k)
          end do
       end do
    end do
  end subroutine calc_3metricdot
  
  subroutine calc_3metricderivdot (gg, dgg, ddgg, kk, dkk, &
       alfa, dalfa, beta, dbeta, ddbeta, dgg_dot)
    REAL*8, intent(in)  :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3)
    REAL*8, intent(in)  :: kk(3,3), dkk(3,3,3)
    REAL*8, intent(in)  :: alfa, dalfa(3)
    REAL*8, intent(in)  :: beta(3), dbeta(3,3), ddbeta(3,3,3)
    REAL*8, intent(out) :: dgg_dot(3,3,3)
    integer :: i,j,k,l
    
    
    
    
    
    
    
    
    do i=1,3
       do j=1,3
          do k=1,3
             dgg_dot(i,j,k) = - 2*dalfa(k) * kk(i,j) - 2*alfa * dkk(i,j,k)
             do l=1,3
                dgg_dot(i,j,k) = dgg_dot(i,j,k) &
                     + dgg(l,j,k) * dbeta(l,i) + gg(l,j) * ddbeta(l,i,k) &
                     + dgg(i,l,k) * dbeta(l,j) + gg(i,l) * ddbeta(l,j,k) &
                     + dbeta(l,k) * dgg(i,j,l) + beta(l) * ddgg(i,j,l,k)
             end do
          end do
       end do
    end do
  end subroutine calc_3metricderivdot
  
  
  
  
  subroutine calc_3metricdot_simple (kk, gg_dot)
    REAL*8, intent(in)  :: kk(3,3)
    REAL*8, intent(out) :: gg_dot(3,3)
    integer :: i,j
    
    
    
    do i=1,3
       do j=1,3
          gg_dot(i,j) = -2 * kk(i,j)
       end do
    end do
  end subroutine calc_3metricdot_simple
  
  subroutine calc_3metricderivdot_simple (dkk, dgg_dot)
    REAL*8, intent(in)  :: dkk(3,3,3)
    REAL*8, intent(out) :: dgg_dot(3,3,3)
    integer :: i,j,k
    
    
    
    do i=1,3
       do j=1,3
          do k=1,3
             dgg_dot(i,j,k) = -2 * dkk(i,j,k)
          end do
       end do
    end do
  end subroutine calc_3metricderivdot_simple
  
  
  







































































  
end module adm_metric
