Skip to content

Commit

Permalink
add misalignment and rotation errors to element 111
Browse files Browse the repository at this point in the history
  • Loading branch information
qianglbl committed Oct 19, 2024
1 parent 7ef1921 commit ad0fba4
Show file tree
Hide file tree
Showing 2 changed files with 151 additions and 1 deletion.
2 changes: 1 addition & 1 deletion src/Appl/BeamLineElem.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1079,7 +1079,7 @@ subroutine getflderrt_BeamLineElem(this,pos,extfld,fldata)
elseif(associated(this%pemfldana)) then
call getfldt_EMfldAna(pos,extfld,this%pemfldana,fldata)
elseif(associated(this%pemfldcart)) then
call getfldt_EMfldCart(pos,extfld,this%pemfldcart,fldata)
call getflderrt_EMfldCart(pos,extfld,this%pemfldcart,fldata)
elseif(associated(this%pemfldcyl)) then
call getfldt_EMfldCyl(pos,extfld,this%pemfldcyl,fldata)
elseif(associated(this%pmult)) then
Expand Down
150 changes: 150 additions & 0 deletions src/Appl/EMfldCart.f90
Original file line number Diff line number Diff line change
Expand Up @@ -1156,4 +1156,154 @@ subroutine getfldt_EMfldCart(pos,extfld,this,fldata)
! call flush_(9)

end subroutine getfldt_EMfldCart

!> get the discrete Ex,Ey,Ez, Bx, By, Bz at given x, y, z, t with errors.
!--------------------------------------------------------------------------------------
subroutine getflderrt_EMfldCart(pos,extfld,this,fldata)
implicit none
include 'mpif.h'
double precision, dimension(4), intent(in) :: pos
type (EMfldCart), intent(in) :: this
type (fielddata), intent(in) :: fldata
double precision, dimension(6), intent(out) :: extfld
double precision :: hx,hy,hz,zz,ab,cd,ef
integer :: ix,jx,kx,ix1,jx1,kx1
double precision :: escale,ww,theta0,tt,tmpcos,tmpsin,z
double complex :: tmptime
real*8, dimension(3) :: poss,tmp,temp
real*8 :: dx,dy,anglex,angley,anglez,zedge

zedge = this%Param(1)
dx = this%Param(7)
dy = this%Param(8)
anglex = this%Param(9)
angley = this%Param(10)
anglez = this%Param(11)
z = pos(3)
temp(1) = pos(1) - dx
temp(2) = pos(2) - dy
tmp(1) = temp(1)*cos(anglez) + temp(2)*sin(anglez)
tmp(2) = -temp(1)*sin(anglez) + temp(2)*cos(anglez)
tmp(3) = pos(3) - zedge
temp(1) = tmp(1)*cos(angley)+tmp(3)*sin(angley)
temp(2) = tmp(2)
temp(3) = -tmp(1)*sin(angley)+tmp(3)*cos(angley)
tmp(1) = temp(1)
tmp(2) = temp(2)*cos(anglex)+temp(3)*sin(anglex)
tmp(3) = -temp(2)*sin(anglex)+temp(3)*cos(anglex)
poss(1:3) = tmp(1:3)
!z = tmp(3)

if((z.gt.this%Param(1)).and.(z.lt.(this%Param(1)+this%Length))) then
escale = this%Param(2)
ww = this%Param(3)*4*asin(1.0)
theta0 = this%Param(4)*asin(1.0)/90
tt = pos(4)
! tmpcos = cos(ww*tt+theta0)*escale
!tmpsin = sin(ww*tt+theta0)*escale
! tmpsin = -sin(ww*tt+theta0)*escale
tmptime = dcmplx(cos(ww*tt+theta0),sin(ww*tt+theta0))*escale
hx = (fldata%XmaxRfgt-fldata%XminRfgt)/fldata%NxIntvRfgt
ix = (poss(1)-fldata%XminRfgt)/hx + 1
ab = ((fldata%XminRfgt-poss(1))+ix*hx)/hx
hy = (fldata%YmaxRfgt-fldata%YminRfgt)/fldata%NyIntvRfgt
jx = (poss(2)-fldata%YminRfgt)/hy + 1
cd = ((fldata%YminRfgt-poss(2))+jx*hy)/hy
hz = (fldata%ZmaxRfgt-fldata%ZminRfgt)/fldata%NzIntvRfgt
zz = poss(3)-this%Param(1)+fldata%ZminRfgt
kx = (zz-fldata%ZminRfgt)/hz + 1
ef = ((fldata%ZminRfgt-zz)+kx*hz)/hz
ix1 = ix + 1
jx1 = jx + 1
kx1 = kx + 1

! print*,"EMCART, ix,jx,kx: ",ix,jx,kx

extfld(1) = dreal((fldata%Exgridt(ix,jx,kx)*ab*cd*ef &
+fldata%Exgridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Exgridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Exgridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Exgridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Exgridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Exgridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Exgridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

extfld(2) = dreal((fldata%Eygridt(ix,jx,kx)*ab*cd*ef &
+fldata%Eygridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Eygridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Eygridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Eygridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Eygridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Eygridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Eygridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

extfld(3) = dreal((fldata%Ezgridt(ix,jx,kx)*ab*cd*ef &
+fldata%Ezgridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Ezgridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Ezgridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Ezgridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Ezgridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Ezgridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Ezgridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

extfld(4) = dreal((fldata%Bxgridt(ix,jx,kx)*ab*cd*ef &
+fldata%Bxgridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Bxgridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Bxgridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Bxgridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Bxgridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Bxgridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Bxgridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

extfld(5) = dreal((fldata%Bygridt(ix,jx,kx)*ab*cd*ef &
+fldata%Bygridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Bygridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Bygridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Bygridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Bygridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Bygridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Bygridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

extfld(6) = dreal((fldata%Bzgridt(ix,jx,kx)*ab*cd*ef &
+fldata%Bzgridt(ix,jx1,kx)*ab*(1.d0-cd)*ef &
+fldata%Bzgridt(ix,jx1,kx1)*ab*(1.d0-cd)*(1.d0-ef) &
+fldata%Bzgridt(ix,jx,kx1)*ab*cd*(1.d0-ef) &
+fldata%Bzgridt(ix1,jx,kx1)*(1.d0-ab)*cd*(1.d0-ef) &
+fldata%Bzgridt(ix1,jx1,kx1)*(1.d0-ab)*(1.d0-cd)*(1.d0-ef)&
+fldata%Bzgridt(ix1,jx1,kx)*(1.d0-ab)*(1.d0-cd)*ef &
+fldata%Bzgridt(ix1,jx,kx)*(1.d0-ab)*cd*ef)*tmptime)

!transform
!for E field
tmp(1) = extfld(1)
tmp(2) = extfld(2)*cos(anglex)-extfld(3)*sin(anglex)
tmp(3) = extfld(2)*sin(anglex)+extfld(3)*cos(anglex)
temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
temp(2) = tmp(2)
temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
extfld(1) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
extfld(2) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
extfld(3) = temp(3)

!for B field
tmp(1) = extfld(4)
tmp(2) = extfld(5)*cos(anglex)-extfld(6)*sin(anglex)
tmp(3) = extfld(5)*sin(anglex)+extfld(6)*cos(anglex)
temp(1) = tmp(1)*cos(angley)-tmp(3)*sin(angley)
temp(2) = tmp(2)
temp(3) = tmp(1)*sin(angley)+tmp(3)*cos(angley)
extfld(4) = temp(1)*cos(anglez) - temp(2)*sin(anglez)
extfld(5) = temp(1)*sin(anglez) + temp(2)*cos(anglez)
extfld(6) = temp(3)

else
extfld = 0.0
endif
! write(9,100)pos(1),pos(2),pos(3),ix*1.0,jx*1.0,kx*1.0,&
! extfld(1:6),fldata%Exgridt(ix,jx,kx),fldata%Eygridt(ix,jx,kx),&
! fldata%Ezgridt(ix,jx,kx)
!100 format(15(1x,e14.6))
! call flush_(9)

end subroutine getflderrt_EMfldCart
end module EMfldCartclass

0 comments on commit ad0fba4

Please sign in to comment.