; FUNCTION createperiodicboundary, x, Box, Boundary, dim = dim
; 
; created by Ken Desmond probably 2011, commented by Eric Weeks Aug 2025
;
; dim = 2 or 3
; 
; Box = size of box [x1,x2,y1,y2] (can also add z1,z2)
; Boundary = mix of rigid 'R' and periodic 'P'
;    for example, Boundary = ['R', 'R'] or ['P', 'R']
;    specifying the X and Y boundary conditions (and Z if needed)
;
; This program takes the original data and adds new data on all
; sides of the original data.  For rigid boundaries, the new data
; is a reflection around the boundary.  For periodic boundaries,
; the new data is the periodic image of the original data.


; first a helper function:

FUNCTION kd_length, a, dim

s = size(a)

IF s[0] EQ 0 THEN return, 0

IF n_elements(dim) EQ 0 THEN BEGIN

   return, s[1]

ENDIF ELSE BEGIN

   IF dim GT s[0] THEN BEGIN
      return, 0
   ENDIF ELSE BEGIN
      return, s[dim]
   ENDELSE

ENDELSE

END
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

FUNCTION createperiodicboundary, x, Box, Boundary, dim = dim



IF ~keyword_set(dim) THEN dim = kd_length(x, 1)
n = kd_length(x, 2)
xp = x


IF dim EQ 2 THEN BEGIN 
   dx = Box[1]-Box[0]
   dy = Box[3]-Box[2]
   
   P = [[-dx, -dy], $
        [ 0 , -dy], $
        [ dx, -dy], $
        [-dx,  0 ], $
        [ dx,  0 ], $
        [-dx,  dy], $
        [ 0 ,  dy], $
        [ dx,  dy]]
   
   
   ;periodic boundary
   FOR i = 0, 7 DO BEGIN
      
      IF (P[0, i] NE 0 AND strcmp(Boundary[0], 'P')) OR (P[0, i] EQ 0) THEN BEGIN
         IF (P[1, i] NE 0 AND strcmp(Boundary[1], 'P')) OR (P[1, i] EQ 0) THEN BEGIN
            
            xp = [[xp],[xp[0,0:n-1]+P[0,i],xp[1,0:n-1]+P[1,i]]]
            
         ENDIF
      ENDIF
      
   ENDFOR
   
   ;If only one boundary is reflexive
   IF strcmp(Boundary[0], 'R')+strcmp(Boundary[1], 'R') EQ 1 THEN BEGIN 
      
      n = kd_length(xp, 2)
      
      IF strcmp(Boundary[0], 'R') THEN BEGIN 
         beta = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(-1.0),xp[1,0:n-1]]]
         beta = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(1.0),xp[1,0:n-1]]]
      ENDIF
      IF strcmp(Boundary[1], 'R') THEN BEGIN 
         beta = xp[1, 0:n-1]-Box[2]
         xp = [[xp],[xp[0,0:n-1],xp[1,0:n-1]+2*beta*(-1.0)]]
         beta = Box[3]-xp[1, 0:n-1]
         xp = [[xp],[xp[0,0:n-1],xp[1,0:n-1]+2*beta*(1.0)]]
      ENDIF 
      
   ENDIF
   
   ;all reflexive boundaries
   IF strcmp(Boundary[0], 'R')+strcmp(Boundary[1], 'R') EQ 2 THEN BEGIN 
      
      ;bottom left
      beta2 = x[1, *]-Box[2]         
      beta1 = x[0, *]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0),xp[1,0:n-1]+2*beta2*(-1.0)]]
      
      ;bottom middle
      beta = x[1, *]-Box[2]
      xp = [[xp],[xp[0,0:n-1],xp[1,0:n-1]+2*beta*(-1.0)]]
      
      ;bottom right
      beta2 = x[1, *]-Box[2]         
      beta1 = Box[1]-x[0, *]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0),xp[1,0:n-1]+2*beta2*(-1.0)]]
      
      ;middle left
      beta = x[0, *]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta*(-1.0),xp[1,0:n-1]]]   
      
      ;middle right 
      beta = Box[1]-x[0, *]
      xp = [[xp],[xp[0,0:n-1]+2*beta*(1.0),xp[1,0:n-1]]]
      
      ;top left 
      beta2 = Box[3]-x[1, *]      
      beta1 = x[0, *]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0),xp[1,0:n-1]+2*beta2*(1.0)]]   
      
      ;top middle
      beta = Box[3]-x[1, *]
      xp = [[xp],[xp[0,0:n-1],xp[1,0:n-1]+2*beta*(1.0)]]
      
      ;top right
      beta2 = Box[3]-x[1, *]      
      beta1 = Box[1]-x[0, *]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0),xp[1,0:n-1]+2*beta2*(1.0)]]
         
   ENDIF

ENDIF

IF dim EQ 3 THEN BEGIN 

   dx = Box[1]-Box[0]
   dy = Box[3]-Box[2]
   dz = Box[5]-Box[4]

   P = make_array(3, 26, value=0.0)

   cnt = 0
   FOR i = -dx, dx, dx DO BEGIN
      FOR j = -dy, dy, dy DO BEGIN
         FOR k = -dz, dz, dz DO BEGIN
            IF abs(i)+abs(j)+abs(k) NE 0 THEN BEGIN
               P[*, cnt] = [i, j, k]
               cnt = cnt+1
            ENDIF
         ENDFOR
      ENDFOR
   ENDFOR


   ;periodic boundary
   FOR i = 0, 25 DO BEGIN
      IF (P[0, i] NE 0 AND strcmp(Boundary[0], 'P')) OR (P[0, i] EQ 0) THEN BEGIN
         IF (P[1, i] NE 0 AND strcmp(Boundary[1], 'P')) OR (P[1, i] EQ 0) THEN BEGIN
            IF (P[2, i] NE 0 AND strcmp(Boundary[2], 'P')) OR (P[2, i] EQ 0) THEN BEGIN
            
               xp = [[xp],[xp[0,0:n-1]+P[0,i], xp[1,0:n-1]+P[1,i], xp[2,0:n-1]+P[2,i]]]

            ENDIF
         ENDIF
      ENDIF      
   ENDFOR


   IF strcmp(Boundary[0], 'R')+strcmp(Boundary[1], 'R')+strcmp(Boundary[2], 'R') EQ 1 THEN BEGIN 

      n = kd_length(xp, 2)
      
      IF strcmp(Boundary[0], 'R') THEN BEGIN 
         beta = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(-1.0), xp[1,0:n-1], xp[2,0:n-1]]]
         beta = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta*( 1.0), xp[1,0:n-1], xp[2,0:n-1]]]
      ENDIF
      IF strcmp(Boundary[1], 'R') THEN BEGIN 
         beta = xp[1, 0:n-1]-Box[2]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*(-1.0), xp[2,0:n-1]]]
         beta = Box[3]-xp[1, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*( 1.0), xp[2,0:n-1]]]
      ENDIF      
      IF strcmp(Boundary[2], 'R') THEN BEGIN 
         beta = xp[2, 0:n-1]-Box[4]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*(-1.0)]]
         beta = Box[5]-xp[2, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*( 1.0)]]
      ENDIF  

   ENDIF

   ;one boundary is confined or periodic
   ;this case is just like 2d all reflexive boundary
   IF strcmp(Boundary[0], 'R')+strcmp(Boundary[1], 'R')+strcmp(Boundary[2], 'R') EQ 2 THEN BEGIN 

      n = kd_length(xp, 2)

      IF strcmp(Boundary[0], 'R') EQ 0 THEN BEGIN

         ;bottom left
         beta2 = xp[1, 0:n-1]-Box[2]         
         beta1 = xp[2, 0:n-1]-Box[4]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta1*(-1.0)]]
      
         ;bottom middle
         beta = xp[1, 0:n-1]-Box[2]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*(-1.0), xp[2,0:n-1]]]
      
         ;bottom right
         beta2 = xp[1, 0:n-1]-Box[2]         
         beta1 = Box[5]-xp[2, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta1*(1.0)]]
      
         ;middle left
         beta = xp[2, 0:n-1]-Box[4]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*(-1.0)]]   
      
         ;middle right 
         beta = Box[5]-xp[2, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*(1.0)]]
      
         ;top left 
         beta2 = Box[3]-xp[1, 0:n-1]      
         beta1 = xp[2, 0:n-1]-Box[4]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta2*(1.0), xp[2,0:n-1]+2*beta1*(-1.0)]]   
      
         ;top middle
         beta = Box[3]-xp[1, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*(1.0), xp[2,0:n-1]]]
      
         ;top right
         beta2 = Box[3]-xp[1, 0:n-1]      
         beta1 = Box[5]-xp[2, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0), xp[1,0:n-1]+2*beta2*(1.0), xp[2,0:n-1]+2*beta1*(1.0)]]


      ENDIF

      IF strcmp(Boundary[1], 'R') EQ 0 THEN BEGIN

         ;bottom left
         beta2 = xp[2, 0:n-1]-Box[4]         
         beta1 = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1], xp[2,0:n-1]+2*beta2*(-1.0)]]
      
         ;bottom middle
         beta = xp[2, 0:n-1]-Box[4]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*(-1.0)]]
      
         ;bottom right
         beta2 = xp[2, 0:n-1]-Box[4]         
         beta1 = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0), xp[1,0:n-1], xp[2,0:n-1]+2*beta2*(-1.0)]]
      
         ;middle left
         beta = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(-1.0), xp[1,0:n-1], xp[2,0:n-1]]]   
      
         ;middle right 
         beta = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(1.0), xp[1,0:n-1], xp[2,0:n-1]]]
      
         ;top left 
         beta2 = Box[5]-xp[2, 0:n-1]      
         beta1 = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1], xp[2,0:n-1]+2*beta2*(1.0)]]   
      
         ;top middle
         beta = Box[5]-xp[2, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1], xp[2,0:n-1]+2*beta*(1.0)]]
      
         ;top right
         beta2 = Box[5]-xp[2, 0:n-1]      
         beta1 = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0), xp[1,0:n-1], xp[2,0:n-1]+2*beta2*(1.0)]]

      ENDIF

      IF strcmp(Boundary[2], 'R') EQ 0 THEN BEGIN

         ;bottom left
         beta2 = xp[1, 0:n-1]-Box[2]         
         beta1 = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]]]
      
         ;bottom middle
         beta = xp[1, 0:n-1]-Box[2]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*(-1.0), xp[2,0:n-1]]]
      
         ;bottom right
         beta2 = xp[1, 0:n-1]-Box[2]         
         beta1 = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]]]
      
         ;middle left
         beta = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(-1.0), xp[1,0:n-1], xp[2,0:n-1]]]   
      
         ;middle right 
         beta = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta*(1.0), xp[1,0:n-1], xp[2,0:n-1]]]
      
         ;top left 
         beta2 = Box[3]-xp[1, 0:n-1]      
         beta1 = xp[0, 0:n-1]-Box[0]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(1.0), xp[2,0:n-1]]]   
      
         ;top middle
         beta = Box[3]-xp[1, 0:n-1]
         xp = [[xp],[xp[0,0:n-1], xp[1,0:n-1]+2*beta*(1.0), xp[2,0:n-1]]]
      
         ;top right
         beta2 = Box[3]-xp[1, 0:n-1]      
         beta1 = Box[1]-xp[0, 0:n-1]
         xp = [[xp],[xp[0,0:n-1]+2*beta1*(1.0), xp[1,0:n-1]+2*beta2*(1.0), xp[2,0:n-1]]]


      ENDIF

   ENDIF

   ;all reflexive boundaries
   IF strcmp(Boundary[0], 'R')+strcmp(Boundary[1], 'R')+strcmp(Boundary[2], 'R') EQ 3 THEN BEGIN 


      ;bottom 
      ;-----------------------------------------------------------------------------------
      ;bottom lower left
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
            
      ;bottom lower middle
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]

      ;bottom lower right
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      
      ;bottom middle left
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = 0    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;bottom middle middle
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = 0    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;bottom middle right
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = 0   
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;bottom upper left
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;bottom upper middle
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
 
      ;bottom upper right
      beta3 = xp[2, 0:n-1]-Box[4]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      ;-----------------------------------------------------------------------------------






      ;center
      ;-----------------------------------------------------------------------------------
      ;center lower left
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center lower middle
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center lower right
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center middle left
      beta3 = 0
      beta2 = 0    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center middle right
      beta3 = 0
      beta2 = 0    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center upper left
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;center upper middle
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      
      ;center upper right
      beta3 = 0
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      ;-----------------------------------------------------------------------------------




      ;top
      ;-----------------------------------------------------------------------------------
      ;top lower left
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top lower middle
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top lower right
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[2]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      
      ;top middle left
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = 0  
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top middle middle
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = 0  
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top middle right
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = 0  
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top upper left
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[0]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   

      ;top upper middle
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = 0
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      
      ;top upper right
      beta3 = xp[2, 0:n-1]-Box[5]
      beta2 = xp[1, 0:n-1]-Box[3]    
      beta1 = xp[0, 0:n-1]-Box[1]
      xp = [[xp],[xp[0,0:n-1]+2*beta1*(-1.0), xp[1,0:n-1]+2*beta2*(-1.0), xp[2,0:n-1]+2*beta2*(-1.0)]]   
      ;-----------------------------------------------------------------------------------


   ENDIF



ENDIF







return, xp

END