From 96ae896517bcb83fa741cd203892cb42a88f0e4f Mon Sep 17 00:00:00 2001
From: Ji Qiang <jqiang@lbl.gov>
Date: Sat, 27 Jul 2024 15:40:29 -0700
Subject: [PATCH] fixed the defect with the new quadratric 3D field
 interpolation

---
 src/Appl/BeamBunch.f90 | 378 ++++++++++++++++++++++++++++-------------
 1 file changed, 264 insertions(+), 114 deletions(-)

diff --git a/src/Appl/BeamBunch.f90 b/src/Appl/BeamBunch.f90
index 2f2b023..d70a000 100644
--- a/src/Appl/BeamBunch.f90
+++ b/src/Appl/BeamBunch.f90
@@ -2795,30 +2795,55 @@ subroutine scatter2_BeamBunch(innp,innx,inny,innz,rays,exg,&
 !                efx = (XminRfg-xx+ix*hxx)*hxxi
 !                efy = (YminRfg-yy+iy*hyy)*hyyi
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) = dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&
@@ -2888,30 +2913,55 @@ subroutine scatter2_BeamBunch(innp,innx,inny,innz,rays,exg,&
                 iy = (yy-YminRfg)*hyyi + 1
                 iy1 = iy+1
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) =extfld(1)+dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&
@@ -3200,30 +3250,55 @@ subroutine scatter20_BeamBunch(innp,rays,tg,gam,chge,mass,&
                 iy = (yy-YminRfg)*hyyi + 1
                 iy1 = iy+1
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) = dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&
@@ -3294,30 +3369,55 @@ subroutine scatter20_BeamBunch(innp,rays,tg,gam,chge,mass,&
                 iy = (yy-YminRfg)*hyyi + 1
                 iy1 = iy+1
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) =extfld(1)+dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&
@@ -5565,30 +5665,55 @@ subroutine scatter2wake_BeamBunch(innp,innx,inny,innz,rays,exg,&
 !                efx = (XminRfg-xx+ix*hxx)*hxxi
 !                efy = (YminRfg-yy+iy*hyy)*hyyi
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) = dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&
@@ -5658,30 +5783,55 @@ subroutine scatter2wake_BeamBunch(innp,innx,inny,innz,rays,exg,&
                 iy = (yy-YminRfg)*hyyi + 1
                 iy1 = iy+1
 
-                ab=(rays(1,n)-XminRfg-(ix-1)*hxx)*hxxi
-                if(ab.le.0.5d0) then
-                   ix2 = ix - 1
+                ab=(xx-XminRfg-(ix-1)*hxx)*hxxi
+                if((ix.eq.1).and.(ab.le.0.5d0)) then
+                   ix2 = ix 
                    wix = 0.75-ab*ab
                    wix1 = (0.5+ab)**2/2
-                   wix2 = (0.5-ab)**2/2
-                else
-                   ix2 = ix + 2
-                   wix = (1.5d0-ab)**2/2
-                   wix1 = 0.75 - (1-ab)*(1-ab)
-                   wix2 = (ab-0.5)**2/2
+                   wix2 = 0.0d0
+                else if((ix.eq.NxIntvRfg).and.(ab.gt.0.5d0)) then
+                    ix2 = ix
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = 0.0d0
+                else 
+                  if(ab.le.0.5d0) then
+                    ix2 = ix - 1
+                    wix = 0.75-ab*ab
+                    wix1 = (0.5+ab)**2/2
+                    wix2 = (0.5-ab)**2/2
+                  else
+                    ix2 = ix + 2
+                    wix = (1.5d0-ab)**2/2
+                    wix1 = 0.75 - (1-ab)*(1-ab)
+                    wix2 = (ab-0.5)**2/2
+                  endif
                 endif
 
-                cd=(rays(3,n)-YminRfg-(iy-1)*hyy)*hyyi
-                if(cd.le.0.5d0) then
-                   iy2 = iy - 1
-                   wiy = 0.75-cd*cd
-                   wiy1 = (0.5+cd)**2/2
-                   wiy2 = (0.5-cd)**2/2
+                cd=(yy-YminRfg-(iy-1)*hyy)*hyyi
+
+                if((iy.eq.1).and.(cd.le.0.5d0)) then
+                    iy2 = iy 
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = 0.0d0
+                else if((iy.eq.NyIntvRfg).and.(cd.gt.0.5d0)) then
+                    iy2 = iy 
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = 0.0d0
                 else
-                   iy2 = iy + 2
-                   wiy = (1.5d0-cd)**2/2
-                   wiy1 = 0.75 - (1-cd)*(1-cd)
-                   wiy2 = (cd-0.5)**2/2
+                  if(cd.le.0.5d0) then
+                    iy2 = iy - 1
+                    wiy = 0.75-cd*cd
+                    wiy1 = (0.5+cd)**2/2
+                    wiy2 = (0.5-cd)**2/2
+                  else
+                    iy2 = iy + 2
+                    wiy = (1.5d0-cd)**2/2
+                    wiy1 = 0.75 - (1-cd)*(1-cd)
+                    wiy2 = (cd-0.5)**2/2
+                  endif
                 endif
 
                 extfld(1) =extfld(1)+dreal((extfld6xyz(1,ix,iy)*wix*wiy+extfld6xyz(1,ix,iy1)*wix*wiy1+&