#include "cctk.h" #include "cctk_Parameters.h" #include "cctk_Arguments.h" #include "cctk_DefineThorn.h" subroutine ADMadjust(CCTK_ARGUMENTS) USE ADM_Scalars implicit none DECLARE_CCTK_ARGUMENTS DECLARE_CCTK_PARAMETERS DECLARE_CCTK_FUNCTIONS integer :: i,j,k,nx,ny,nz,kk,ll,adjlevel,flag_adj_guts CCTK_REAL :: zero,o13 CCTK_REAL :: alp_ijk CCTK_REAL :: gxx_ijk,gxy_ijk,gxz_ijk,gyy_ijk,gyz_ijk,gzz_ijk CCTK_REAL :: kxx_ijk,kxy_ijk,kxz_ijk,kyy_ijk,kyz_ijk,kzz_ijk CCTK_REAL :: gxx_ijk2,gxy_ijk2,gxz_ijk2,gyy_ijk2,gyz_ijk2,gzz_ijk2 CCTK_REAL :: agxx,agxy,agxz,agyy,agyz,agzz CCTK_REAL :: akxx,akxy,akxz,akyy,akyz,akzz CCTK_REAL :: aPxx,aPxy,aPxz,aPyy,aPyz,aPzz CCTK_REAL :: aRxx,aRxy,aRxz,aRyy,aRyz,aRzz CCTK_REAL :: aQxx(3),aQxy(3),aQxz(3),aQyy(3),aQyz(3),aQzz(3) CCTK_REAL :: aSxx(3),aSxy(3),aSxz(3),aSyy(3),aSyz(3),aSzz(3) CCTK_REAL :: appxx(3),appxy(3),appxz(3),appyy(3),appyz(3),appzz(3) CCTK_REAL :: arrxx(3),arrxy(3),arrxz(3),arryy(3),arryz(3),arrzz(3) CCTK_REAL :: aqqxx(3,3),aqqxy(3,3),aqqxz(3,3) & ,aqqyy(3,3),aqqyz(3,3),aqqzz(3,3) CCTK_REAL :: assxx(3,3),assxy(3,3),assxz(3,3) & ,assyy(3,3),assyz(3,3),asszz(3,3) CCTK_REAL :: deltaf CCTK_REAL :: trK,ugij(3,3),alpdx(3),ahamdx(3),amomdx(3,3) #include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing_declare.h" c Declarations for macros from ADM Evolution (Hisaaki Shinkai) #include "macro_adjustedADM/DH_declare.h" #include "macro_adjustedADM/DMOM_declare.h" #include "macro_adjustedADM/CDMOM_declare.h" c Declarations for macros from CactusEinstein #include "EinsteinBase/ADMMacros/src/macro/DA_declare.h" #include "EinsteinBase/ADMMacros/src/macro/CHR2_declare.h" c for the adjusted coef. ... c flag_adj_guts=0 (use previous step: same with constr), c =1 (use current step) flag_adj_guts=1 if (CCTK_EQUALS(adjust_type,'standard')) then continue else zero = 0.0D0 o13=1.0D0/3.0D0 nx = cctk_lsh(1) ny = cctk_lsh(2) nz = cctk_lsh(3) #include "EinsteinBase/ADMMacros/src/macro/ADM_Spacing.h" adjlevel=1 c zero out ahamdx = zero amomdx = zero agxx=zero;agxy=zero;agxz=zero;agyy=zero;agyz=zero;agzz=zero akxx=zero;akxy=zero;akxz=zero;akyy=zero;akyz=zero;akzz=zero aPxx=zero;aPxy=zero;aPxz=zero;aPyy=zero;aPyz=zero;aPzz=zero aQxx=zero;aQxy=zero;aQxz=zero;aQyy=zero;aQyz=zero;aQzz=zero aRxx=zero;aRxy=zero;aRxz=zero;aRyy=zero;aRyz=zero;aRzz=zero aSxx=zero;aSxy=zero;aSxz=zero;aSyy=zero;aSyz=zero;aSzz=zero appxx=zero;appxy=zero;appxz=zero appyy=zero;appyz=zero;appzz=zero aqqxx=zero;aqqxy=zero;aqqxz=zero aqqyy=zero;aqqyz=zero;aqqzz=zero arrxx=zero;arrxy=zero;arrxz=zero arryy=zero;arryz=zero;arrzz=zero assxx=zero;assxy=zero;assxz=zero assyy=zero;assyz=zero;asszz=zero c----------------------------------------------------------------------- c prepare constraints -> removed in 2002-June version c----------------------------------------------------------------------- c call adjADMConstraints(CCTK-ARGUMENTS) c----------------------------------------------------------------------- c adjustments c----------------------------------------------------------------------- do k=2,nz-1 do j=2,ny-1 do i=2,nx-1 if(flag_adj_guts.eq.0) then alp_ijk=ADM_alp_p(i,j,k) gxx_ijk=ADM_gxx_p(i,j,k) gxy_ijk=ADM_gxy_p(i,j,k) gxz_ijk=ADM_gxz_p(i,j,k) gyy_ijk=ADM_gyy_p(i,j,k) gyz_ijk=ADM_gyz_p(i,j,k) gzz_ijk=ADM_gzz_p(i,j,k) kxx_ijk=ADM_kxx_p(i,j,k) kxy_ijk=ADM_kxy_p(i,j,k) kxz_ijk=ADM_kxz_p(i,j,k) kyy_ijk=ADM_kyy_p(i,j,k) kyz_ijk=ADM_kyz_p(i,j,k) kzz_ijk=ADM_kzz_p(i,j,k) else alp_ijk=alp(i,j,k) gxx_ijk=gxx(i,j,k) gxy_ijk=gxy(i,j,k) gxz_ijk=gxz(i,j,k) gyy_ijk=gyy(i,j,k) gyz_ijk=gyz(i,j,k) gzz_ijk=gzz(i,j,k) kxx_ijk=kxx(i,j,k) kxy_ijk=kxy(i,j,k) kxz_ijk=kxz(i,j,k) kyy_ijk=kyy(i,j,k) kyz_ijk=kyz(i,j,k) kzz_ijk=kzz(i,j,k) endif if (CCTK_EQUALS(adjust_type,'original')) then c== originalADM type aRxx=-0.25D0*alp_ijk*gxx_ijk aRxy=-0.25D0*alp_ijk*gxy_ijk aRxz=-0.25D0*alp_ijk*gxz_ijk aRyy=-0.25D0*alp_ijk*gyy_ijk aRyz=-0.25D0*alp_ijk*gyz_ijk aRzz=-0.25D0*alp_ijk*gzz_ijk else if (CCTK_EQUALS(adjust_type,'Detweiler')) then adjlevel=2 c== Detweiler type if(i.eq.1.and.j.eq.1.and.k.eq.1) & call CCTK_INFO("adjusted ADM: adjust_type=Detweiler") #include "macro_adjustedADM/detweiler_P.h" #include "macro_adjustedADM/detweiler_S.h" #include "macro_adjustedADM/detweiler_ss.h" c== Detweiler R trK = ugij(1,1)*Kxx_ijk + ugij(2,2)*Kyy_ijk + ugij(3,3)*Kzz_ijk & + 2.0D0*(ugij(1,2)*Kxy_ijk + ugij(1,3)*Kxz_ijk + ugij(2,3)*Kyz_ijk) aRxx= adjust_kappa*alp_ijk**3*(kxx_ijk-o13*trK*kxx_ijk) aRxy= adjust_kappa*alp_ijk**3*(kxy_ijk-o13*trK*kxy_ijk) aRxz= adjust_kappa*alp_ijk**3*(kxz_ijk-o13*trK*kxz_ijk) aRyy= adjust_kappa*alp_ijk**3*(kyy_ijk-o13*trK*kyy_ijk) aRyz= adjust_kappa*alp_ijk**3*(kyz_ijk-o13*trK*kyz_ijk) aRzz= adjust_kappa*alp_ijk**3*(kzz_ijk-o13*trK*kzz_ijk) else if((CCTK_EQUALS(adjust_type,'DetweilerP')).or. & (CCTK_EQUALS(adjust_type,'P1'))) then c== Detweiler P-type #include "macro_adjustedADM/detweiler_P.h" else if (CCTK_EQUALS(adjust_type,'DetweilerS')) then c== Detweiler S-type #include "macro_adjustedADM/detweiler_S.h" else if((CCTK_EQUALS(adjust_type,'Detweilerss')).or. & (CCTK_EQUALS(adjust_type,'s_1'))) then c== Detweiler ss-type adjlevel=2 #include "macro_adjustedADM/detweiler_ss.h" else if((CCTK_EQUALS(adjust_type,'simpleDet')).or. & (CCTK_EQUALS(adjust_type,'P2'))) then c== simplified Detweiler aPxx=-adjust_kappa*alp_ijk*gxx_ijk aPxy=-adjust_kappa*alp_ijk*gxy_ijk aPxz=-adjust_kappa*alp_ijk*gxz_ijk aPyy=-adjust_kappa*alp_ijk*gyy_ijk aPyz=-adjust_kappa*alp_ijk*gyz_ijk aPzz=-adjust_kappa*alp_ijk*gzz_ijk else if((CCTK_EQUALS(adjust_type,'constantRxx')).or. & (CCTK_EQUALS(adjust_type,'R2'))) then c== constantRxx aRxx=-adjust_kappa*alp_ijk aRyy=-adjust_kappa*alp_ijk aRzz=-adjust_kappa*alp_ijk else if (CCTK_EQUALS(adjust_type,'P3')) then c== P-3 aPxx=-adjust_kappa aPyy=-adjust_kappa aPzz=-adjust_kappa else if (CCTK_EQUALS(adjust_type,'P4')) then c== P-4 aPxx=-adjust_kappa*gxx_ijk aPxy=-adjust_kappa*gxy_ijk aPxz=-adjust_kappa*gxz_ijk aPyy=-adjust_kappa*gyy_ijk aPyz=-adjust_kappa*gyz_ijk aPzz=-adjust_kappa*gzz_ijk else if (CCTK_EQUALS(adjust_type,'P5')) then c== P-5 aPxx=-adjust_kappa*gxx_ijk aPyy=-adjust_kappa*gyy_ijk aPzz=-adjust_kappa*gzz_ijk else if (CCTK_EQUALS(adjust_type,'Q1')) then c== Q-1 do kk=1,3 aQxx(kk)=adjust_kappa*alp_ijk*gxx_ijk enddo else if (CCTK_EQUALS(adjust_type,'Q2')) then c== Q-2 else if (CCTK_EQUALS(adjust_type,'Q3')) then c== Q-3 else if (CCTK_EQUALS(adjust_type,'Q4')) then c== Q-4 else if (CCTK_EQUALS(adjust_type,'R1')) then c== R-1 else if (CCTK_EQUALS(adjust_type,'R3')) then c== R-3 else if (CCTK_EQUALS(adjust_type,'S1')) then c== S-1 else if (CCTK_EQUALS(adjust_type,'S2')) then c== S-2 else if (CCTK_EQUALS(adjust_type,'p_1')) then c== p-1 adjlevel=2 else if (CCTK_EQUALS(adjust_type,'p_2')) then c== p-2 adjlevel=2 else if (CCTK_EQUALS(adjust_type,'p_3')) then c== p-3 adjlevel=2 else if (CCTK_EQUALS(adjust_type,'q_1')) then c== q-1 adjlevel=2 do ll=1,3 do kk=1,3 if(kk.eq.ll) then aqqxx(kk,ll)=adjust_kappa*alp_ijk*gxx_ijk aqqxy(kk,ll)=adjust_kappa*alp_ijk*gxy_ijk aqqxz(kk,ll)=adjust_kappa*alp_ijk*gxz_ijk aqqyy(kk,ll)=adjust_kappa*alp_ijk*gyy_ijk aqqyz(kk,ll)=adjust_kappa*alp_ijk*gyz_ijk aqqzz(kk,ll)=adjust_kappa*alp_ijk*gzz_ijk endif enddo enddo else if (CCTK_EQUALS(adjust_type,'q_2')) then c== q-2 adjlevel=2 do ll=1,3 do kk=1,3 if(kk.eq.ll) then aqqxx(kk,ll)=-adjust_kappa*alp_ijk*gxx_ijk aqqyy(kk,ll)=-adjust_kappa*alp_ijk*gyy_ijk aqqzz(kk,ll)=-adjust_kappa*alp_ijk*gzz_ijk endif enddo enddo else if (CCTK_EQUALS(adjust_type,'r_1')) then c== r-1 adjlevel=2 do kk=1,3 c if(kk.eq.1) then arrxx(kk)= adjust_kappa*alp_ijk*gxx_ijk arrxy(kk)= adjust_kappa*alp_ijk*gxy_ijk arrxz(kk)= adjust_kappa*alp_ijk*gxz_ijk arryy(kk)= adjust_kappa*alp_ijk*gyy_ijk arryz(kk)= adjust_kappa*alp_ijk*gyz_ijk arrzz(kk)= adjust_kappa*alp_ijk*gzz_ijk c endif enddo else if (CCTK_EQUALS(adjust_type,'r_2')) then c== r-2 adjlevel=2 do kk=1,3 c if(kk.eq.1) then arrxx(kk)=-adjust_kappa*alp_ijk arryy(kk)=-adjust_kappa*alp_ijk arrzz(kk)=-adjust_kappa*alp_ijk c endif enddo else if (CCTK_EQUALS(adjust_type,'r_3')) then c== r-3 adjlevel=2 do kk=1,3 c if(kk.eq.1) then arrxx(kk)= adjust_kappa*alp_ijk*gxx_ijk arrxy(kk)= adjust_kappa*alp_ijk*gxy_ijk arrxz(kk)= adjust_kappa*alp_ijk*gxz_ijk arryy(kk)= adjust_kappa*alp_ijk*gyy_ijk arryz(kk)= adjust_kappa*alp_ijk*gyz_ijk arrzz(kk)= adjust_kappa*alp_ijk*gzz_ijk c endif enddo else if (CCTK_EQUALS(adjust_type,'s_2')) then c== s-2 adjlevel=2 do kk=1,3 if(kk.eq.ll) then assxx(kk,ll)=-adjust_kappa*alp_ijk*gxx_ijk assyy(kk,ll)=-adjust_kappa*alp_ijk*gyy_ijk asszz(kk,ll)=-adjust_kappa*alp_ijk*gzz_ijk endif enddo else if (CCTK_EQUALS(adjust_type,'s_3')) then c== s-3 do ll=1,3 do kk=1,3 if(kk.eq.ll) then assxx(kk,ll)=-adjust_kappa*alp_ijk*gxx_ijk assyy(kk,ll)=-adjust_kappa*alp_ijk*gyy_ijk asszz(kk,ll)=-adjust_kappa*alp_ijk*gzz_ijk endif enddo enddo adjlevel=2 else call CCTK_WARN(0,"Specify adjust_type: adjustedADM.F") endif c----------------------------------------------------------------------- c delivatives of constraints if(adjlevel.eq.2) then #include "macro_adjustedADM/DH_guts.h" ahamdx(1)=DH_DXDH ahamdx(2)=DH_DYDH ahamdx(3)=DH_DZDH c temporal stock in order to use macro/CHR2_guts.h gxx_ijk2=gxx(i,j,k) gxy_ijk2=gxy(i,j,k) gxz_ijk2=gxz(i,j,k) gyy_ijk2=gyy(i,j,k) gyz_ijk2=gyz(i,j,k) gzz_ijk2=gzz(i,j,k) gxx(i,j,k)=gxx_ijk gxy(i,j,k)=gxy_ijk gxz(i,j,k)=gxz_ijk gyy(i,j,k)=gyy_ijk gyz(i,j,k)=gyz_ijk gzz(i,j,k)=gzz_ijk #include "EinsteinBase/ADMMacros/src/macro/CHR2_guts.h" gxx(i,j,k)=gxx_ijk2 gxy(i,j,k)=gxy_ijk2 gxz(i,j,k)=gxz_ijk2 gyy(i,j,k)=gyy_ijk2 gyz(i,j,k)=gyz_ijk2 gzz(i,j,k)=gzz_ijk2 #include "macro_adjustedADM/DMOM_guts.h" #include "macro_adjustedADM/CDMOM_guts.h" amomdx(1,1)=CD_DXDMOMX amomdx(1,2)=CD_DXDMOMY amomdx(1,3)=CD_DXDMOMZ amomdx(2,1)=CD_DYDMOMX amomdx(2,2)=CD_DYDMOMY amomdx(2,3)=CD_DYDMOMZ amomdx(3,1)=CD_DZDMOMX amomdx(3,2)=CD_DZDMOMY amomdx(3,3)=CD_DZDMOMZ endif #include "macro_adjustedADM/DH_undefine.h" #include "macro_adjustedADM/DMOM_undefine.h" #include "macro_adjustedADM/CDMOM_undefine.h" #include "EinsteinBase/ADMMacros/src/macro/CHR2_undefine.h" c----------------------------------------------------------------------- c adjustment in gij equations agxx=aPxx*aham(i,j,k)+aQxx(1)*amomx(i,j,k) & +aQxx(2)*amomy(i,j,k) & +aQxx(3)*amomz(i,j,k) agxy=aPxy*aham(i,j,k)+aQxy(1)*amomx(i,j,k) & +aQxy(2)*amomy(i,j,k) & +aQxy(3)*amomz(i,j,k) agxz=aPxz*aham(i,j,k)+aQxz(1)*amomx(i,j,k) & +aQxz(2)*amomy(i,j,k) & +aQxz(3)*amomz(i,j,k) agyy=aPyy*aham(i,j,k)+aQyy(1)*amomx(i,j,k) & +aQyy(2)*amomy(i,j,k) & +aQyy(3)*amomz(i,j,k) agyz=aPyz*aham(i,j,k)+aQyz(1)*amomx(i,j,k) & +aQyz(2)*amomy(i,j,k) & +aQyz(3)*amomz(i,j,k) agzz=aPzz*aham(i,j,k)+aQzz(1)*amomx(i,j,k) & +aQzz(2)*amomy(i,j,k) & +aQzz(3)*amomz(i,j,k) if(adjlevel.eq.2) then do kk=1,3 agxx=agxx+appxx(kk)*ahamdx(kk) agxy=agxy+appxy(kk)*ahamdx(kk) agxz=agxz+appxz(kk)*ahamdx(kk) agyy=agyy+appyy(kk)*ahamdx(kk) agyz=agyz+appyz(kk)*ahamdx(kk) agzz=agzz+appzz(kk)*ahamdx(kk) end do do kk=1,3 do ll=1,3 agxx=agxx+aqqxx(kk,ll)*amomdx(kk,ll) agxy=agxy+aqqxy(kk,ll)*amomdx(kk,ll) agxz=agxz+aqqxz(kk,ll)*amomdx(kk,ll) agyy=agyy+aqqyy(kk,ll)*amomdx(kk,ll) agyz=agyz+aqqyz(kk,ll)*amomdx(kk,ll) agzz=agzz+aqqzz(kk,ll)*amomdx(kk,ll) end do end do endif adms_gxx(i,j,k) = adms_gxx(i,j,k) + agxx adms_gxy(i,j,k) = adms_gxy(i,j,k) + agxy adms_gxz(i,j,k) = adms_gxz(i,j,k) + agxz adms_gyy(i,j,k) = adms_gyy(i,j,k) + agyy adms_gyz(i,j,k) = adms_gyz(i,j,k) + agyz adms_gzz(i,j,k) = adms_gzz(i,j,k) + agzz c adjustment in Kij equations akxx=aRxx*aham(i,j,k)+aSxx(1)*amomx(i,j,k) & +aSxx(2)*amomy(i,j,k) & +aSxx(3)*amomz(i,j,k) akxy=aRxy*aham(i,j,k)+aSxy(1)*amomx(i,j,k) & +aSxy(2)*amomy(i,j,k) & +aSxy(3)*amomz(i,j,k) akxz=aRxz*aham(i,j,k)+aSxz(1)*amomx(i,j,k) & +aSxz(2)*amomy(i,j,k) & +aSxz(3)*amomz(i,j,k) akyy=aRyy*aham(i,j,k)+aSyy(1)*amomx(i,j,k) & +aSyy(2)*amomy(i,j,k) & +aSyy(3)*amomz(i,j,k) akyz=aRyz*aham(i,j,k)+aSyz(1)*amomx(i,j,k) & +aSyz(2)*amomy(i,j,k) & +aSyz(3)*amomz(i,j,k) akzz=aRzz*aham(i,j,k)+aSzz(1)*amomx(i,j,k) & +aSzz(2)*amomy(i,j,k) & +aSzz(3)*amomz(i,j,k) if(adjlevel.eq.2) then do kk=1,3 akxx=akxx+arrxx(kk)*ahamdx(kk) akxy=akxy+arrxy(kk)*ahamdx(kk) akxz=akxz+arrxz(kk)*ahamdx(kk) akyy=akyy+arryy(kk)*ahamdx(kk) akyz=akyz+arryz(kk)*ahamdx(kk) akzz=akzz+arrzz(kk)*ahamdx(kk) end do do kk=1,3 do ll=1,3 akxx=akxx+assxx(kk,ll)*amomdx(kk,ll) akxy=akxy+assxy(kk,ll)*amomdx(kk,ll) akxz=akxz+assxz(kk,ll)*amomdx(kk,ll) akyy=akyy+assyy(kk,ll)*amomdx(kk,ll) akyz=akyz+assyz(kk,ll)*amomdx(kk,ll) akzz=akzz+asszz(kk,ll)*amomdx(kk,ll) end do end do endif adms_kxx(i,j,k) = adms_kxx(i,j,k) + akxx adms_kxy(i,j,k) = adms_kxy(i,j,k) + akxy adms_kxz(i,j,k) = adms_kxz(i,j,k) + akxz adms_kyy(i,j,k) = adms_kyy(i,j,k) + akyy adms_kyz(i,j,k) = adms_kyz(i,j,k) + akyz adms_kzz(i,j,k) = adms_kzz(i,j,k) + akzz end do end do end do endif return end subroutine ADMadjust c####################################################################### c----------------------------------------------------------------------- function deltaf(i,j) implicit none integer :: i,j CCTK_REAL :: deltaf if(i.eq.j) then deltaf=1.0D0 else deltaf=0.0D0 endif end function deltaf c####################################################################### c----------------------------------------------------------------------- c B_{ij}=inverse of A_{ij} c input .... Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz c output ... Bxx,Bxy,Bxz,Byx,Byy,Byz,Bzx,Bzy,Bzz subroutine inverse33real(Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz, & Bxx,Bxy,Bxz,Byx,Byy,Byz,Bzx,Bzy,Bzz) implicit none c input-output CCTK_REAL :: Axx,Axy,Axz,Ayx,Ayy,Ayz,Azx,Azy,Azz, & Bxx,Bxy,Bxz,Byx,Byy,Byz,Bzx,Bzy,Bzz c locally defined CCTK_REAL :: bunbo c bunbo= b11*b22*b33 - b11*b23*b32 - b21*b12*b33 c & + b21*b13*b32 + b31*b12*b23 - b31*b13*b22 bunbo= Axx*Ayy*Azz - Axx*Ayz*Azy - Ayx*Axy*Azz & + Ayx*Axz*Azy + Azx*Axy*Ayz - Azx*Axz*Ayy Bxx= (Ayy*Azz - Ayz*Azy)/bunbo Bxy=-(Axy*Azz - Axz*Azy)/bunbo Bxz= (Axy*Ayz - Axz*Ayy)/bunbo Byx=-(Ayx*Azz - Ayz*Azx)/bunbo Byy= (Axx*Azz - Axz*Azx)/bunbo Byz=-(Axx*Ayz - Axz*Ayx)/bunbo Bzx= (Ayx*Azy - Ayy*Azx)/bunbo Bzy=-(Axx*Azy - Axy*Azx)/bunbo Bzz= (Axx*Ayy - Axy*Ayx)/bunbo c$$$[m22 m33 - m23 m32 m12 m33 - m13 m32 m12 m23 - m13 m22] c$$$[----------------- , - ----------------- , -----------------] c$$$[ %1 %1 %1 ] c$$$[ ] c$$$[ m21 m33 - m23 m31 m11 m33 - m13 m31 m11 m23 - m13 m21] c$$$[- ----------------- , ----------------- , - -----------------] c$$$[ %1 %1 %1 ] c$$$[ ] c$$$[m21 m32 - m22 m31 m11 m32 - m12 m31 m11 m22 - m12 m21] c$$$[----------------- , - ----------------- , -----------------] c$$$[ %1 %1 %1 ] c$$$ c$$$%1 := m11 m22 m33 - m11 m23 m32 - m21 m12 m33 c$$$ c$$$ + m21 m13 m32 + m31 m12 m23 - m31 m13 m22 return end c#######################################################################