Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fixes to exact drift when pt is large and beta << 1 #1077

Merged
merged 2 commits into from Feb 25, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
77 changes: 39 additions & 38 deletions src/twiss.f90
Expand Up @@ -7139,7 +7139,7 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te)
use twissbeamfi, only : beta, gamma, dtbyds
use twisslfi
use matrices, only : EYE
use math_constfi, only : zero, two, three
use math_constfi, only : zero, one, two, three
implicit none
!----------------------------------------------------------------------*
! Purpose: *
Expand Down Expand Up @@ -7172,54 +7172,55 @@ SUBROUTINE tmdrf(fsec,ftrk,orbit,fmap,dl,ek,re,te)
py = orbit(4)
pt = orbit(6)

csq = 1 + 2*pt*beta + pt**2 - px**2 - py**2
csq = 1 + 2*pt/beta + pt**2 - px**2 - py**2
l_pz = dl / sqrt(csq)
c3sq = csq**(3d0/2d0)
c52sq =csq**(5d0/2d0)
re(1,2) = dl/sqrt(csq) + dl*px**2/c3sq;
re(1,4) = dl*px*py/c3sq;
re(1,6) = dl*px*(-beta - pt)/c3sq;
re(3,2) = dl*px*py/c3sq;
re(3,4) = dl/sqrt(csq) + dl*py**2/c3sq;
re(3,6) = dl*py*(-beta - pt)/c3sq;
re(5,2) = -dl*px*(beta + pt)/c3sq;
re(5,4) = -dl*py*(beta + pt)/c3sq;
re(5,6) = -dl/sqrt(csq) - dl*(-beta - pt)*(beta + pt)/c3sq;

re(1,2) = dl/sqrt(csq) + dl*px**2/c3sq
re(1,4) = dl*px*py/c3sq
re(1,6) = dl*px*(-pt - 1d0/beta)/c3sq
re(3,2) = dl*px*py/c3sq
re(3,4) = dl/sqrt(csq) + dl*py**2/c3sq
re(3,6) = dl*py*(-pt - 1d0/beta)/c3sq
re(5,2) = -dl*px*(beta + pt)/c3sq
re(5,4) = -dl*py*(beta + pt)/c3sq
re(5,6) = -dl/sqrt(csq) - dl*(beta + pt)*(-pt - 1d0/beta)/c3sq

if (fsec) then
te(1,2,2) = 3d0*dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**3d0/(2d0*c52sq)
te(1,2,4) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(1,2,6) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + dl*px**2d0*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
te(1,4,2) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(1,4,4) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(1,4,6) = dl*px*py*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
te(1,6,2) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*(-beta - pt)/(2d0*c52sq)
te(1,6,4) = 3d0*dl*px*py*(-beta - pt)/(2d0*c52sq)
te(1,6,6) = -dl*px/(2d0*csq**(3d0/2d0)) + dl*px*(-3d0*beta - 3d0*pt)*(-beta - pt)/(2d0*c52sq)
te(3,2,2) = dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(3,2,4) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(3,2,6) = dl*px*py*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
te(3,4,2) = dl*px/(2d0*csq**(3d0/2d0)) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(3,4,4) = 3d0*dl*py/(2d0*csq**(3d0/2d0)) + 3d0*dl*py**3d0/(2d0*c52sq)
te(3,4,6) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + dl*py**2d0*(-3d0*beta - 3d0*pt)/(2d0*c52sq)
te(3,6,2) = 3d0*dl*px*py*(-beta - pt)/(2d0*c52sq)
te(3,6,4) = dl*(-beta - pt)/(2d0*csq**(3d0/2d0)) + 3d0*dl*py**2d0*(-beta - pt)/(2d0*c52sq)
te(3,6,6) = -dl*py/(2d0*csq**(3d0/2d0)) + dl*py*(-3d0*beta - 3d0*pt)*(-beta - pt)/(2d0*c52sq)
te(5,2,2) = -dl*(beta + pt)/(2d0*csq**(3d0/2d0)) - 3d0*dl*px**2d0*(beta + pt)/(2d0*c52sq)
te(1,2,2) = 3d0*dl*px/(2d0*c3sq) + 3d0*dl*px**3d0/(2d0*c52sq)
te(1,2,4) = dl*py/(2d0*c3sq) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(1,2,6) = dl*(-pt - 1d0/beta)/(2d0*c3sq) + dl*px**2d0*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(1,4,2) = dl*py/(2d0*c3sq) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(1,4,4) = dl*px/(2d0*c3sq) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(1,4,6) = dl*px*py*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(1,6,2) = dl*(-pt - 1d0/beta)/(2d0*c3sq) + 3d0*dl*px**2d0*(-pt - 1d0/beta)/(2d0*c52sq)
te(1,6,4) = 3d0*dl*px*py*(-pt - 1d0/beta)/(2d0*c52sq)
te(1,6,6) = -dl*px/(2d0*c3sq) + dl*px*(-3d0*pt - 3d0/beta)*(-pt - 1d0/beta)/(2d0*c52sq)
te(3,2,2) = dl*py/(2d0*c3sq) + 3d0*dl*px**2d0*py/(2d0*c52sq)
te(3,2,4) = dl*px/(2d0*c3sq) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(3,2,6) = dl*px*py*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(3,4,2) = dl*px/(2d0*c3sq) + 3d0*dl*px*py**2d0/(2d0*c52sq)
te(3,4,4) = 3d0*dl*py/(2d0*c3sq) + 3d0*dl*py**3d0/(2d0*c52sq)
te(3,4,6) = dl*(-pt - 1d0/beta)/(2d0*c3sq) + dl*py**2d0*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(3,6,2) = 3d0*dl*px*py*(-pt - 1d0/beta)/(2d0*c52sq)
te(3,6,4) = dl*(-pt - 1d0/beta)/(2d0*c3sq) + 3d0*dl*py**2d0*(-pt - 1d0/beta)/(2d0*c52sq)
te(3,6,6) = -dl*py/(2d0*c3sq) + dl*py*(-3d0*pt - 3d0/beta)*(-pt - 1d0/beta)/(2d0*c52sq)
te(5,2,2) = -dl*(beta + pt)/(2d0*c3sq) - 3d0*dl*px**2d0*(beta + pt)/(2d0*c52sq)
te(5,2,4) = -3d0*dl*px*py*(beta + pt)/(2d0*c52sq)
te(5,2,6) = -dl*px/(2d0*csq**(3d0/2d0)) - dl*px*(-3d0*beta - 3d0*pt)*(beta + pt)/(2d0*c52sq)
te(5,2,6) = -dl*px/(2d0*c3sq) - dl*px*(beta + pt)*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(5,4,2) = -3d0*dl*px*py*(beta + pt)/(2d0*c52sq)
te(5,4,4) = -dl*(beta + pt)/(2d0*csq**(3d0/2d0)) - 3d0*dl*py**2d0*(beta + pt)/(2d0*c52sq)
te(5,4,6) = -dl*py/(2d0*csq**(3d0/2d0)) - dl*py*(-3d0*beta - 3d0*pt)*(beta + pt)/(2d0*c52sq)
te(5,6,2) = -dl*px/(2d0*csq**(3d0/2d0)) - 3d0*dl*px*(-beta - pt)*(beta + pt)/(2d0*c52sq)
te(5,6,4) = -dl*py/(2d0*csq**(3d0/2d0)) - 3d0*dl*py*(-beta - pt)*(beta + pt)/(2d0*c52sq)
te(5,6,6) = -dl*(-beta - pt)/csq**(3d0/2d0) + dl*(beta + pt)/(2d0*csq**(3d0/2d0)) &
- dl*(-3d0*beta - 3d0*pt)*(-beta - pt)*(beta + pt)/(2d0*c52sq)
te(5,4,4) = -dl*(beta + pt)/(2d0*c3sq) - 3d0*dl*py**2d0*(beta + pt)/(2d0*c52sq)
te(5,4,6) = -dl*py/(2d0*c3sq) - dl*py*(beta + pt)*(-3d0*pt - 3d0/beta)/(2d0*c52sq)
te(5,6,2) = -dl*px/(2d0*c3sq) - 3d0*dl*px*(beta + pt)*(-pt - 1d0/beta)/(2d0*c52sq)
te(5,6,4) = -dl*py/(2d0*c3sq) - 3d0*dl*py*(beta + pt)*(-pt - 1d0/beta)/(2d0*c52sq)
te(5,6,6) = dl*(beta + pt)/(2d0*c3sq) - dl*(-pt - 1d0/beta)/c3sq - &
dl*(beta + pt)*(-3d0*pt - 3d0/beta)*(-pt - 1d0/beta)/(2d0*c52sq)
endif

orbit(1) = orbit(1) + px*l_pz
orbit(3) = orbit(3) + py*l_pz
orbit(5) = orbit(5) + (dl*beta - (beta + pt) * l_pz)
orbit(5) = orbit(5) + (dl/beta - (1d0/beta + pt) * l_pz)
else

re(1,2) = dl
Expand Down
52 changes: 26 additions & 26 deletions tests/test-twiss-exact/madx_exact.tfs.ref
Expand Up @@ -4,9 +4,9 @@
@ PARTICLE %06s "PROTON"
@ MASS %le 0.93827208816
@ CHARGE %le 1.00000000000
@ ENERGY %le 2.00000000000
@ PC %le 1.76625181913
@ GAMMA %le 2.13157784958
@ ENERGY %le 1.00000000000
@ PC %le 0.34589808988
@ GAMMA %le 1.06578892479
@ KBUNCH %le 1.00000000000
@ BCURRENT %le 0.00000000000
@ SIGE %le 0.00100000000
Expand All @@ -24,16 +24,16 @@
@ Q2 %le 0.21860189884
@ DQ1 %le 0.00000000000
@ DQ2 %le 0.00000000000
@ DXMAX %le 0.00883657376
@ DYMAX %le 0.17673147524
@ DXMAX %le 0.02892764750
@ DYMAX %le 0.57855295009
@ XCOMAX %le 0.01000200560
@ YCOMAX %le 0.20004011206
@ BETXMAX %le 101.04031624715
@ BETYMAX %le 52.06009815221
@ XCORMS %le 0.00638635935
@ YCORMS %le 0.12772718707
@ DXRMS %le 0.00564222194
@ DYRMS %le 0.11284443887
@ DXRMS %le 0.01847053077
@ DYRMS %le 0.36941061531
@ DELTAP %le 0.00000000000
@ SYNCH_1 %le 0.00000000000
@ SYNCH_2 %le 0.00000000000
Expand All @@ -42,24 +42,24 @@
@ SYNCH_5 %le 0.00000000000
@ SYNCH_6 %le 0.00000000000
@ SYNCH_8 %le 0.00000000000
@ DQMIN %le 0.00000158802
@ DQMIN %le 0.00000000000
@ DQMIN_PHASE %le 0.00000004999
@ TITLE %08s "no-title"
@ ORIGIN %16s "5.07.00 Linux 64"
@ DATE %08s "25/08/21"
@ TIME %08s "00.30.01"
* NAME S X Y PT BETX BETY
$ %s %le %le %le %le %le %le
"SEQ$START" 0.00000000000 0.00000000000 0.00000000000 0.00000000000 1.00000000000 2.00000000000
"DRIFT_0" 1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
"M1" 1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
"DRIFT_1" 3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
"M2" 3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
"DRIFT_2" 5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
"M3" 5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
"DRIFT_3" 7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
"M4" 7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
"DRIFT_4" 9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
"M5" 9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
"DRIFT_5" 10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
"SEQ$END" 10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
@ ORIGIN %16s "5.08.00 Linux 64"
@ DATE %08s "02/02/22"
@ TIME %08s "14.40.18"
* NAME S X Y T PT BETX BETY
$ %s %le %le %le %le %le %le %le
"SEQ$START" 0.00000000000 0.00000000000 0.00000000000 0.00000000000 0.00000000000 1.00000000000 2.00000000000
"DRIFT_0" 1.00000000000 0.00100020056 0.02000401121 -0.00057982489 0.00000000000 2.00040316247 2.50060098152
"M1" 1.00000000000 0.00100020056 0.02000401121 -0.00057982489 0.00000000000 2.00040316247 2.50060098152
"DRIFT_1" 3.00000000000 0.00300060168 0.06001203362 -0.00173947466 0.00000000000 10.00362846224 6.50540883370
"M2" 3.00000000000 0.00300060168 0.06001203362 -0.00173947466 0.00000000000 10.00362846224 6.50540883370
"DRIFT_2" 5.00000000000 0.00500100280 0.10002005603 -0.00289912443 0.00000000000 26.01007906179 14.51502453805
"M3" 5.00000000000 0.00500100280 0.10002005603 -0.00289912443 0.00000000000 26.01007906179 14.51502453805
"DRIFT_3" 7.00000000000 0.00700140392 0.14002807844 -0.00405877420 0.00000000000 50.01975496110 26.52944809458
"M4" 7.00000000000 0.00700140392 0.14002807844 -0.00405877420 0.00000000000 50.01975496110 26.52944809458
"DRIFT_4" 9.00000000000 0.00900180504 0.18003610086 -0.00521842397 0.00000000000 82.03265616019 42.54867950329
"M5" 9.00000000000 0.00900180504 0.18003610086 -0.00521842397 0.00000000000 82.03265616019 42.54867950329
"DRIFT_5" 10.00000000000 0.01000200560 0.20004011206 -0.00579824886 0.00000000000 101.04031624715 52.06009815221
"SEQ$END" 10.00000000000 0.01000200560 0.20004011206 -0.00579824886 0.00000000000 101.04031624715 52.06009815221
60 changes: 30 additions & 30 deletions tests/test-twiss-exact/ptc.tfs.ref
Expand Up @@ -4,9 +4,9 @@
@ PARTICLE %06s "PROTON"
@ MASS %le 0.93827208816
@ CHARGE %le 1.00000000000
@ ENERGY %le 2.00000000000
@ PC %le 1.76625181913
@ GAMMA %le 2.13157784958
@ ENERGY %le 1.00000000000
@ PC %le 0.34589808988
@ GAMMA %le 1.06578892479
@ KBUNCH %le 1.00000000000
@ BCURRENT %le 0.00000000000
@ SIGE %le 0.00100000000
Expand Down Expand Up @@ -37,24 +37,24 @@
@ BETA12MIN %le 0.00000000000
@ BETA12MAX %le 0.00000002002
@ BETA13MIN %le 0.00000000000
@ BETA13MAX %le 0.00002567481
@ BETA13MAX %le 0.00016736176
@ BETA21MIN %le 0.00000000000
@ BETA21MAX %le 0.00000004005
@ BETA22MIN %le 2.00000000000
@ BETA22MAX %le 52.06009815221
@ BETA23MIN %le 0.00000000000
@ BETA23MAX %le 0.01026992554
@ BETA23MAX %le 0.06694470321
@ BETA31MIN %le 0.00000000000
@ BETA31MAX %le 0.00012837407
@ BETA31MAX %le 0.00083680879
@ BETA32MIN %le 0.00000000000
@ BETA32MAX %le 0.02567481386
@ BETA32MAX %le 0.16736175803
@ BETA33MIN %le 5.00000000000
@ BETA33MAX %le 6.59915529106
@ DISP1MIN %le -0.00506703206
@ BETA33MAX %le 1089.23251036158
@ DISP1MIN %le -0.01293683725
@ DISP1MAX %le 0.00000000000
@ DISP2MIN %le 0.00000000000
@ DISP2MAX %le 0.00000000000
@ DISP3MIN %le -0.10134064113
@ DISP3MIN %le -0.25873674500
@ DISP3MAX %le 0.00000000000
@ DISP4MIN %le 0.00000000000
@ DISP4MAX %le 0.00000000000
Expand All @@ -68,7 +68,7 @@
@ PXCORMS %le 0.00100000000
@ YCORMS %le 0.12772718707
@ PYCORMS %le 0.02000000000
@ TCORMS %le 0.00145006893
@ TCORMS %le 0.00370222756
@ PTCORMS %le 0.00000000000
@ XCOMIN %le 0.00000000000
@ XCOMAX %le 0.01000200560
Expand All @@ -78,26 +78,26 @@
@ YCOMAX %le 0.20004011206
@ PYCOMIN %le 0.02000000000
@ PYCOMAX %le 0.02000000000
@ TCOMIN %le -0.00227102748
@ TCOMIN %le -0.00579824886
@ TCOMAX %le -0.00000000000
@ PTCOMIN %le 0.00000000000
@ PTCOMAX %le 0.00000000000
@ TITLE %08s "no-title"
@ ORIGIN %16s "5.07.00 Linux 64"
@ DATE %08s "25/08/21"
@ TIME %08s "00.30.01"
* NAME S X Y PT BETX BETY
$ %s %le %le %le %le %le %le
"SEQ$START" 0.00000000000 0.00000000000 0.00000000000 0.00000000000 1.00000000000 2.00000000000
"DRIFT_0" 1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
"M1" 1.00000000000 0.00100020056 0.02000401121 0.00000000000 2.00040316247 2.50060098152
"DRIFT_1" 3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
"M2" 3.00000000000 0.00300060168 0.06001203362 0.00000000000 10.00362846224 6.50540883370
"DRIFT_2" 5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
"M3" 5.00000000000 0.00500100280 0.10002005603 0.00000000000 26.01007906179 14.51502453805
"DRIFT_3" 7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
"M4" 7.00000000000 0.00700140392 0.14002807844 0.00000000000 50.01975496110 26.52944809458
"DRIFT_4" 9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
"M5" 9.00000000000 0.00900180504 0.18003610086 0.00000000000 82.03265616019 42.54867950329
"DRIFT_5" 10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
"SEQ$END" 10.00000000000 0.01000200560 0.20004011206 0.00000000000 101.04031624715 52.06009815221
@ ORIGIN %16s "5.08.00 Linux 64"
@ DATE %08s "02/02/22"
@ TIME %08s "14.42.25"
* NAME S X Y T PT BETX BETY
$ %s %le %le %le %le %le %le %le
"SEQ$START" 0.00000000000 0.00000000000 0.00000000000 -0.00000000000 0.00000000000 1.00000000000 2.00000000000
"DRIFT_0" 1.00000000000 0.00100020056 0.02000401121 -0.00057982489 0.00000000000 2.00040316247 2.50060098152
"M1" 1.00000000000 0.00100020056 0.02000401121 -0.00057982489 0.00000000000 2.00040316247 2.50060098152
"DRIFT_1" 3.00000000000 0.00300060168 0.06001203362 -0.00173947466 0.00000000000 10.00362846224 6.50540883370
"M2" 3.00000000000 0.00300060168 0.06001203362 -0.00173947466 0.00000000000 10.00362846224 6.50540883370
"DRIFT_2" 5.00000000000 0.00500100280 0.10002005603 -0.00289912443 0.00000000000 26.01007906179 14.51502453805
"M3" 5.00000000000 0.00500100280 0.10002005603 -0.00289912443 0.00000000000 26.01007906179 14.51502453805
"DRIFT_3" 7.00000000000 0.00700140392 0.14002807844 -0.00405877420 0.00000000000 50.01975496110 26.52944809458
"M4" 7.00000000000 0.00700140392 0.14002807844 -0.00405877420 0.00000000000 50.01975496110 26.52944809458
"DRIFT_4" 9.00000000000 0.00900180504 0.18003610086 -0.00521842397 0.00000000000 82.03265616019 42.54867950329
"M5" 9.00000000000 0.00900180504 0.18003610086 -0.00521842397 0.00000000000 82.03265616019 42.54867950329
"DRIFT_5" 10.00000000000 0.01000200560 0.20004011206 -0.00579824886 0.00000000000 101.04031624715 52.06009815221
"SEQ$END" 10.00000000000 0.01000200560 0.20004011206 -0.00579824886 0.00000000000 101.04031624715 52.06009815221
4 changes: 2 additions & 2 deletions tests/test-twiss-exact/test-twiss-exact.madx
Expand Up @@ -7,14 +7,14 @@ m4:marker, at=7;
m5:marker, at=9;
endsequence;
set, format="13.11f";
beam,particle=proton,energy=2;
beam,particle=proton,energy=1;
use, sequence=seq;
option, sympl=false;
px0 = 0.001;
py0 = 0.02;
betx0 = 1;
bety0 = 2;
select, flag = twiss, column=name,s,x,y,pt,betx,bety;
select, flag = twiss, column=name,s,x,y,t, pt,betx,bety;

twiss, file="madx.tfs", betx=betx0, bety=bety0, px=px0, py=py0, pt=pt0;
twiss, file="madx_exact.tfs", betx=betx0, bety=bety0, px=px0, py=py0, exact;
Expand Down