Routines for setting up evolution with MoL, for initializing the level set

function and evolving it. Have to add re-parametrization of the level
set function before it can be used.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/EHFinder/trunk@3 2a26948c-0e4f-0410-aee8-f1d3e353619c

7 files changed, 520 insertions, 0 deletions

diff --git a/interface.ccl b/interface.ccl
index 6efef5f..ece47c1 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -1,2 +1,25 @@
 # Interface definition for thorn EHFinder
 # $Header$
+
+implements: ehfinder
+inherits: grid einstein
+
+USES INCLUDE: Einstein.h
+USES INCLUDE: MoL.h
+
+public:
+
+CCTK_REAL level_set type=GF TimeLevels=2
+{
+  f
+} "The scalar level set function that defines the null surface"
+
+CCTK_REAL slevel_set type=GF TimeLevels=1
+{
+  sf
+} "Source for the level set function"
+
+CCTK_REAL dlevel_set type=GF TimeLevels=1
+{
+  dfx, dfy, dfz
+}
diff --git a/param.ccl b/param.ccl
index 004bbdf..9425b58 100644
--- a/param.ccl
+++ b/param.ccl
@@ -1,2 +1,50 @@
 # Parameter definitions for thorn EHFinder
 # $Header$
+
+private:
+
+KEYWORD initial_f "Initial surface choice"
+{
+  "sphere"     :: "spherical surface"
+  "ellipsoid"  :: "ellipsoidal surface"
+  "cassini"    :: "ovaloid of cassini"
+} "sphere"
+
+REAL initial_rad "Initial radius of surface"
+{
+(0.0:            :: "Positive please"
+} 1.0
+
+REAL initial_a "Initial a coefficient of ellipsoid"
+{
+(0.0:            :: "Positive please"
+} 1.0
+
+REAL initial_b "Initial b coefficient of ellipsoid"
+{
+(0.0:            :: "Positive please"
+} 1.0
+
+REAL initial_c "Initial c coefficient of ellipsoid"
+{
+(0.0:            :: "Positive please"
+} 1.0
+
+real cas_a "Initial a coefficient of ovaloid of cassini"
+{
+:                :: "Any number (negative and positive are equivalent)"
+} 2.0
+
+real cas_b "Initial b coefficient of ovaloid of cassini"
+{
+:                :: "Any number (negative and positive are equivalent)"
+} 2.05
+
+
+BOOLEAN normalize "normalize the derivatives of the level set function"
+{
+} "no"
+
+BOOLEAN one_sided "Use one sided differences everywhere"
+{
+} "no"
diff --git a/schedule.ccl b/schedule.ccl
index d3b44f6..39eae35 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -1,2 +1,28 @@
 # Schedule definitions for thorn EHFinder
 # $Header$
+
+STORAGE: level_set
+STORAGE: slevel_set
+STORAGE: dlevel_set
+
+schedule EHFinder_Init at CCTK_POSTINITIAL
+{
+  LANG: Fortran
+  SYNC: level_set
+} "Setup the initial surface"
+
+schedule EHFinder_MoLRegister in MoL_Register
+{
+  LANG: Fortran
+} "Register evolution variables"
+
+schedule EHFinder_Sources in MoL_CalcRHS
+{
+  LANG: Fortran
+} "Calculate the source terms"
+
+schedule GROUP EHFinder_PostStep in MoL_PostStep
+{
+  LANG: Fortran
+  SYNC: level_set
+} "Schedule syncing of level set"
diff --git a/src/EHFinder_Init.F90 b/src/EHFinder_Init.F90
new file mode 100644
index 0000000..406a1ee
--- /dev/null
+++ b/src/EHFinder_Init.F90
@@ -0,0 +1,42 @@
+! Initialisation of the level set function
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Init(CCTK_ARGUMENTS)
+
+  implicit none
+
+  integer :: ierr
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+!  print*,cctk_dim,cctk_lsh,cctk_gsh
+!  print*,Xlevel_set0,Xlevel_set1,Xlevel_set2
+   
+  ierr = CCTK_Equals( initial_f, 'sphere' )
+  if (ierr .eq. 1) then
+!    f = sqrt(x**2+y**2+z**2) - initial_rad
+!    f = x**2+y**2+z**2 - initial_rad**2
+    f = sqrt(x**2+y**2+z**2)**3 - initial_rad**3
+  end if
+  ierr = CCTK_Equals( initial_f, 'ellipsoid' )
+  if (ierr .eq. 1) then
+    f = sqrt(x**2/initial_a**2+y**2/initial_b**2+z**2/initial_c**2) - 1.0
+  end if
+  ierr = CCTK_Equals( initial_f, 'cassini' )
+  if (ierr .eq. 1) then
+    f = (x**2+y**2+z**2)**2 + cas_a**4 - &
+        2*cas_a**2*(x**2 - (y**2+z**2)) - cas_b**4
+  end if
+!  print*,f(:,1,1)
+!  print*,f(:,1,1)
+!  print*,initial_rad
+
+  return
+end subroutine EHFinder_Init
+
+
diff --git a/src/EHFinder_Sources.F90 b/src/EHFinder_Sources.F90
new file mode 100644
index 0000000..48e946a
--- /dev/null
+++ b/src/EHFinder_Sources.F90
@@ -0,0 +1,311 @@
+! Calculation of the sources for the level set function.
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_Sources(CCTK_ARGUMENTS)
+
+  implicit none
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  CCTK_INT :: i, j, k, nx, ny, nz, ngx, ngy, ngz, ierr
+  CCTK_INT :: ix1, ix2, iy1, iy2, iz1, iz2
+  CCTK_INT :: jx1, jx2, jy1, jy2, jz1, jz2
+!  CCTK_REAL, allocatable, dimension(:,:,:) :: dfx, dfy, dfz, dfsq
+  CCTK_REAL, allocatable, dimension(:,:,:) :: dfsq
+  CCTK_REAL :: idx, idy, idz
+  CCTK_REAL :: a, b, c
+  CCTK_REAL :: g4tt, g4tx, g4ty, g4tz, g4xx, g4xy, g4xz, g4yy, g4yz, g4zz
+  CCTK_REAL :: idetg, alp2, tmp1, tmp2, tmp3
+  CCTK_REAL :: dfsq, ratio
+
+  nx = cctk_lsh(1)
+  ny = cctk_lsh(2)
+  nz = cctk_lsh(3)
+  
+  ngx = cctk_nghostzones(1)
+  ngy = cctk_nghostzones(2)
+  ngz = cctk_nghostzones(3)
+
+!  print*,nx,ny,nz
+!  print*,cctk_delta_space
+
+! allocate storage for derivatives of the level set function
+
+!  allocate ( dfx(nx,ny,nz), dfy(nx,ny,nz), dfz(nx,ny,nz), dfsq(nx,ny,nz), stat=ierr )
+  allocate ( dfsq(nx,ny,nz), stat=ierr )
+  if ( ierr .gt. 0 ) then
+    call CCTK_WARN(0,'Unable to allocate work arrays')
+  end if
+
+! calculate 1/(2*delta) in each direction
+
+  idx = 1.0 / ( 2.0 * cctk_delta_space(1) )
+  idy = 1.0 / ( 2.0 * cctk_delta_space(2) )
+  idz = 1.0 / ( 2.0 * cctk_delta_space(3) )
+
+  if ( cctk_bbox(1) .eq. 1 ) then
+    ix1 = 0
+    jx1 = 1
+  else
+    ix1 = ngx
+    jx1 = ngx
+  end if
+  if ( cctk_bbox(2) .eq. 1 ) then
+    ix2 = 0
+    jx2 = 1
+  else
+    ix2 = ngx
+    jx2 = ngx
+  end if
+
+  if ( cctk_bbox(3) .eq. 1 ) then
+    iy1 = 0
+    jy1 = 1
+  else 
+    iy1 = ngy
+    jy1 = ngy
+  endif
+  if ( cctk_bbox(4) .eq. 1 ) then
+    iy2 = 0
+    jy2 = 1
+  else
+    iy2 = ngy
+    jy2 = ngy
+  endif
+
+  if ( cctk_bbox(5) .eq. 1 ) then
+    iz1 = 0
+    jz1 = 1
+  else 
+    iz1 = ngz
+    jz1 = ngz
+  endif
+  if ( cctk_bbox(6) .eq. 1 ) then
+    iz2 = 0
+    jz2 = 1
+  else
+    iz2 = ngz
+    jz2 = ngz
+  endif
+
+! calculate one sided derivatives for the outer boundary points
+! for x
+
+  if ( cctk_bbox(1) .eq. 1 ) then
+    dfx(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
+                                  -3.0*f(1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
+                                   4.0*f(2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+                                       f(3, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+  endif
+  if ( cctk_bbox(2) .eq. 1 ) then
+    dfx(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = idx * ( &
+                                   3.0*f(nx, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+                                   4.0*f(nx-1, 1+iy1:ny-iy2, 1+iz1:nz-iz2) + &
+                                       f(nx-2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+  endif
+
+! for y
+
+  if ( cctk_bbox(3) .eq. 1 ) then
+    dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) = idy * ( &
+                                  -3.0*f(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2) + &
+                                   4.0*f(1+ix1:nx-ix2, 2, 1+iz1:nz-iz2) - &
+                                       f(1+ix1:nx-ix2, 3, 1+iz1:nz-iz2) )
+  endif
+  if ( cctk_bbox(4) .eq. 1 ) then
+    dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) = idy * ( &
+                                   3.0*f(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2) - &
+                                   4.0*f(1+ix1:nx-ix2, ny-1, 1+iz1:nz-iz2) + &
+                                       f(1+ix1:nx-ix2, ny-2, 1+iz1:nz-iz2) )
+  endif
+
+! for z
+
+  if ( cctk_bbox(5) .eq. 1 ) then
+    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) = idz * ( &
+                                  -3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1) + &
+                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2) - &
+                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 3) )
+  endif
+  if ( cctk_bbox(6) .eq. 1 ) then
+    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) = idz * ( &
+                                   3.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz) - &
+                                   4.0*f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-1) + &
+                                       f(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz-2) )
+  endif
+
+! calculate normal centered derivatives for interior points
+
+  if ( one_sided ) then
+    do k = 1+iz1, nz-iz2
+      do j = 1+iy1, ny-iy2
+        do i = 1+jx1, nx-jx2
+          if ( f(i,j,k) - f(i-1,j,k) .ne. 0 ) then
+            ratio = abs ( f(i+1,j,k) - f(i,j,k) ) / ( f(i,j,k) - f(i-1,j,k) )
+            if ( ratio .lt. 1 ) then
+              dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
+            else
+              dfx(i,j,k) = idx * ( 3.0*f(i,j,k) - 4.0*f(i-1,j,k) + f(i-2,j,k) )
+            end if
+          else
+            dfx(i,j,k) = idx * ( -3.0*f(i,j,k) + 4.0*f(i+1,j,k) - f(i+2,j,k) )
+          end if
+        end do
+      end do
+    end do
+          
+    do k = 1+iz1, nz-iz2
+      do j = 1+jy1, ny-jy2
+        do i = 1+ix1, nx-ix2
+          if ( f(i,j,k) - f(i,j-1,k) .ne. 0 ) then
+            ratio = abs ( f(i,j+1,k) - f(i,j,k) ) / ( f(i,j,k) - f(i,j-1,k) )
+            if ( ratio .lt. 1 ) then
+              dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
+            else
+              dfy(i,j,k) = idy * ( 3.0*f(i,j,k) - 4.0*f(i,j-1,k) + f(i,j-2,k) )
+            end if
+          else
+            dfy(i,j,k) = idy * ( -3.0*f(i,j,k) + 4.0*f(i,j+1,k) - f(i,j+2,k) )
+          end if
+        end do
+      end do
+    end do
+            
+    do k = 1+jz1, nz-jz2
+      do j = 1+iy1, ny-iy2
+        do i = 1+ix1, nx-ix2
+          if ( f(i,j,k) - f(i,j,k-1) .ne. 0 ) then
+            ratio = abs ( f(i,j,k+1) - f(i,j,k) ) / ( f(i,j,k) - f(i,j,k-1) )
+            if ( ratio .lt. 1 ) then
+              dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
+            else
+              dfz(i,j,k) = idz * ( 3.0*f(i,j,k) - 4.0*f(i,j,k-1) + f(i,j,k-2) )
+            end if
+          else
+            dfz(i,j,k) = idz * ( -3.0*f(i,j,k) + 4.0*f(i,j,k+1) - f(i,j,k+2) )
+          end if
+        end do
+      end do
+    end do
+
+  else
+
+    dfx(1+jx1:nx-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) = &
+                  idx * ( f(2+jx1:nx+1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) - &
+                          f(jx1:nx-1-jx2, 1+iy1:ny-iy2, 1+iz1:nz-iz2) )
+    dfy(1+ix1:nx-ix2, 1+jy1:ny-jy2, 1+iz1:nz-iz2) = &
+                  idy * ( f(1+ix1:nx-ix2, 2+jy1:ny+1-jy2, 1+iz1:nz-iz2) - &
+                          f(1+ix1:nx-ix2, jy1:ny-1-jy2, 1+iz1:nz-iz2) )
+    dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1+jz1:nz-jz2) = &
+                  idz * ( f(1+ix1:nx-ix2, 1+iy1:ny-iy2, 2+jz1:nz+1-jz2) - &
+                          f(1+ix1:nx-ix2, 1+iy1:ny-iy2, jz1:nz-1-jz2) )
+  end if
+
+!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
+!    print*,dfx(1,1+iy1:ny-iy2, 1+iz1:nz-iz2)
+!    print*,dfx(nx,1+iy1:ny-iy2, 1+iz1:nz-iz2)
+!    print*
+!    print*,dfy(1+ix1:nx-ix2, 1, 1+iz1:nz-iz2)
+!    print*,dfy(1+ix1:nx-ix2, ny, 1+iz1:nz-iz2)
+!    print*
+!    print*,dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, 1)
+!    print*,dfz(1+ix1:nx-ix2, 1+iy1:ny-iy2, nz)
+!    print*
+!  end if
+! note I multiply by the square of the lapse in the components of
+! the 4 metric. This does not effect the solution to the quadratic
+! equation, but avoids potential problems of dividing with the square 
+! of the lapse in regions where it is close to zero.
+
+        
+  if ( normalize ) then
+    dfsq = 1.0 / sqrt(dfx**2+dfy**2+dfz**2)
+    dfx = dfx*dfsq
+    dfy = dfy*dfsq
+    dfz = dfz*dfsq
+  end if
+
+  do k = 1+iz1, nz-iz2
+    do j = 1+iy1, ny-iy2
+      do i = 1+ix1, nx-ix2
+        alp2 = alp(i,j,k)**2
+        g4tt = -1.0
+        g4tx = betax(i,j,k)
+        g4ty = betay(i,j,k)
+        g4tz = betaz(i,j,k)
+
+        tmp1 = gyy(i,j,k)*gzz(i,j,k) - gyz(i,j,k)**2
+        tmp2 = gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k)
+        tmp3 = gxy(i,j,k)*gyz(i,j,k) - gxz(i,j,k)*gyy(i,j,k)
+
+        idetg = 1.0 / ( gxx(i,j,k)*tmp1 + gxy(i,j,k)*tmp2 + gxz(i,j,k)*tmp3 )
+
+        g4xx = tmp1 * idetg
+        g4xy = tmp2 * idetg
+        g4xz = tmp3 * idetg
+        g4yy = ( gxx(i,j,k)*gzz(i,j,k) - gxz(i,j,k)**2 ) * idetg
+        g4yz = ( gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k) ) * idetg
+        g4zz = ( gxx(i,j,k)*gyy(i,j,k) - gxy(i,j,k)**2 ) * idetg
+        
+        g4xx = alp2*g4xx - betax(i,j,k)**2
+        g4xy = alp2*g4xy - betax(i,j,k)*betay(i,j,k)
+        g4xz = alp2*g4xz - betax(i,j,k)*betaz(i,j,k)
+        g4yy = alp2*g4yy - betay(i,j,k)**2
+        g4yz = alp2*g4yz - betay(i,j,k)*betaz(i,j,k)
+        g4zz = alp2*g4zz - betaz(i,j,k)**2
+      
+        a = g4tt
+        b = 2.0 * ( g4tx * dfx(i,j,k) + &
+                    g4ty * dfy(i,j,k) + &
+                    g4tz * dfz(i,j,k) )
+        c = g4xx * dfx(i,j,k)**2 + &
+            g4yy * dfy(i,j,k)**2 + &
+            g4zz * dfz(i,j,k)**2 + &
+            2.0 * ( g4xy * dfx(i,j,k) * dfy(i,j,k) + &
+                    g4xz * dfx(i,j,k) * dfz(i,j,k) + &
+                    g4yz * dfy(i,j,k) * dfz(i,j,k) )
+
+        if ( b .lt. 0. ) then
+          sf(i,j,k) = - ( -b + sqrt(b**2 - 4*a*c) ) / ( 2.0*a )
+        else
+          sf(i,j,k) = - 2.0*c / ( -b - sqrt(b**2 - 4*a*c) )
+        end if
+!        if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
+!          if ( abs(b) .gt. 0.00000001 ) then
+!            print*,i,j,k
+!            print*,dfx(i,j,k),dfy(i,j,k),dfz(i,j,k)
+!            print*,g4tx,g4ty,g4tz
+!            print*,idetg
+!            print*,g4xx,g4yy,g4zz
+!            print*,g4xy,g4xz,g4yz
+!            print*,a,b,c
+!            print*,sf(i,j,k)
+!            pause
+!          endif
+!        end if
+      end do
+    end do
+  end do 
+!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
+!    print*,sf
+!    print*
+!  end if
+!  print*,cctk_bbox
+!  print*,cctk_nghostzones
+
+!  deallocate ( dfx, dfy, dfz, dfsq )
+  deallocate ( dfsq )
+!  if ( CCTK_MyProc(cctkGH) .eq. 0 ) then
+!    call CCTK_WARN(0,'stopping')
+!  end if
+  
+  return
+end subroutine EHFinder_Sources
+
+
diff --git a/src/MoL_Init.F90 b/src/MoL_Init.F90
new file mode 100644
index 0000000..ed6585c
--- /dev/null
+++ b/src/MoL_Init.F90
@@ -0,0 +1,59 @@
+! Registration of variables with MoL
+! $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+subroutine EHFinder_MoLRegister(CCTK_ARGUMENTS)
+
+  implicit none
+
+  CCTK_INT :: ierr, ierr_cum, varindex, rhsindex
+
+  DECLARE_CCTK_PARAMETERS
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+
+  ierr_cum = 0
+  call CCTK_VarIndex(varindex, 'ehfinder::f')
+  call CCTK_VarIndex(rhsindex, 'ehfinder::sf')
+  call MoL_RegisterVar(ierr, varindex, rhsindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gxx')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gxy')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gxz')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gyy')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gyz')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::gzz')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::alp')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::betax')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::betay')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  call CCTK_VarIndex(varindex, 'einstein::betaz')
+  call MoL_RegisterPrimitive(ierr, varindex)
+  ierr_cum = ierr_cum + ierr
+  print*,'MoLRegister ', ierr_cum
+  if ( ierr_cum .gt. 0 ) then
+    call CCTK_WARN(0,'Problems registering variables with MoL')
+  end if
+end subroutine EHFinder_MoLRegister
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..b237039
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,11 @@
+# Main make.code.defn file for thorn EHFinder
+# $Header$
+
+# Source files in this directory
+SRCS = EHFinder_Init.F90 \
+       MoL_Init.F90 \
+       EHFinder_Sources.F90
+
+# Subdirectories containing source files
+SUBDIRS = 
+