Merge pull request #858 from tpersson/fixThinRF

Fixes a number of reported issues
MethodicalAcceleratorDesign · Nov 15, 2019 · d11e158 · d11e158
2 parents d326d0a + d86065b
commit d11e158
Show file tree

Hide file tree

Showing 4 changed files with 43 additions and 20 deletions.
diff --git a/doc/latexuguide/elements.tex b/doc/latexuguide/elements.tex
@@ -900,6 +900,11 @@ \section{RF Cavity}
     in particular to \texttt{TIME} and \texttt{TOTALPATH} switches,
     concerning RF behaviour of cavities in PTC.
 
+  \item The charge is not taken into account by the RF-cavity.
+
+  \item The RF-cavity has to be thin for TRACKING in MAD-X. This is handled by 
+  \texttt{MAKETHIN} but note that it is always sliced into a single slice.  
+
 \end{itemize}
 
 The \texttt{RFCAVITY} can also have attributes that only become active in

diff --git a/src/mad_mkthin.cpp b/src/mad_mkthin.cpp
@@ -2674,6 +2674,12 @@ sequence* SequenceList::slice_sequence(const std::string slice_style,sequence* t
     if (theSeqElList.current_node() == thick_sequ->end)
     {
       break;
+    } 
+    if(theSeqElList.current_node()->p_elem!=nullptr){
+      if(strcmp(theSeqElList.current_node()->p_elem->base_type->name, "rfcavity")==0 &&
+        find_element(theSeqElList.current_node()->p_elem->name, sliced_seq->cavities) == nullptr){
+        add_to_el_list(&theSeqElList.current_node()->p_elem, 0, sliced_seq->cavities, 0);
+      }
     }
     theSeqElList.current_node(theSeqElList.current_node()->next); // set current_node
   }
@@ -2695,5 +2701,7 @@ sequence* SequenceList::slice_sequence(const std::string slice_style,sequence* t
   put_sequ(thick_sequ); // Slicing done for this sequence. Add to list of sequences sliced
   if(MaTh::Verbose) std::cout << __FILE__ << " " << __FUNCTION__ << " line " << std::setw(4) << __LINE__ << " before print theSeqElList" << std::endl;
   if(MaTh::Verbose) theSeqElList.Print(); // print final list
+
+
   return sliced_seq;
 } // slice_sequence
diff --git a/src/trrun.f90 b/src/trrun.f90
@@ -3560,7 +3560,7 @@ subroutine trsol(track,ktrack,dxt,dyt)
 
            !---- Radiation loss at entrance (step.eq.1) and exit (step.eq.3)
            if ((step.eq.1).or.(step.eq.3)) then
-              if (radiate) then
+              if (radiate .and. elrad .gt. zero) then
                  !---- Full damping.
                  if (damp) then
                     curv = sqrt(dxt(i)**2 + dyt(i)**2) / elrad;
@@ -3630,7 +3630,7 @@ subroutine trsol(track,ktrack,dxt,dyt)
               dyt(i) = pyf_ - track(4,i);
 
               if ((step.eq.1).or.(step.eq.3)) then
-                 if (radiate) then
+                 if (radiate .and. elrad .gt. zero) then
                     !---- Full damping.
                     if (damp) then
                        curv = sqrt(dxt(i)**2 + dyt(i)**2) / length;

diff --git a/src/twiss.f90 b/src/twiss.f90
@@ -7943,15 +7943,20 @@ SUBROUTINE tmsol_th(ftrk,orbit,fmap,ek,re,te)
 
   !---- Half radiation effect at entry.
   if (radiate .and. ftrk) then
-     kx = ((sk**2)*orbit(1)-sk*orbit(4))*elrad;
-     ky = ((sk**2)*orbit(3)+sk*orbit(2))*elrad;
-     rfac = (arad * gamma**3 / three) * (kx**2 + ky**2) / elrad;
-     pt = orbit(6);
-     bet_sqr = (pt*pt + two*pt/bet0 + one) / (one/bet0 + pt)**2;
-     f_damp_t = sqrt(one + rfac*(rfac - two) / bet_sqr);
-     orbit(2) = orbit(2) * f_damp_t;
-     orbit(4) = orbit(4) * f_damp_t;
-     orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
+    if(elrad .eq. zero) then 
+      call fort_warn('TWCPGO: ','Radiation effects ignored for solenoid '// &
+                         'with l=0, lrad=0 and radiate=true')
+    else
+       kx = ((sk**2)*orbit(1)-sk*orbit(4))*elrad;
+       ky = ((sk**2)*orbit(3)+sk*orbit(2))*elrad;
+       rfac = (arad * gamma**3 / three) * (kx**2 + ky**2) / elrad;
+       pt = orbit(6);
+       bet_sqr = (pt*pt + two*pt/bet0 + one) / (one/bet0 + pt)**2;
+       f_damp_t = sqrt(one + rfac*(rfac - two) / bet_sqr);
+       orbit(2) = orbit(2) * f_damp_t;
+       orbit(4) = orbit(4) * f_damp_t;
+       orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
+    endif
   endif
 
   !---- First-order terms.
@@ -7977,15 +7982,20 @@ SUBROUTINE tmsol_th(ftrk,orbit,fmap,ek,re,te)
 
   !---- Half radiation effect at exit.
   if (radiate .and. ftrk) then
-     kx = ((sk**2)*orbit(1)-sk*orbit(4))*elrad;
-     ky = ((sk**2)*orbit(3)+sk*orbit(2))*elrad;
-     rfac = (arad * gamma**3 / three) * (kx**2 + ky**2) / elrad;
-     pt = orbit(6);
-     bet_sqr = (pt*pt + two*pt/bet0 + one) / (one/bet0 + pt)**2;
-     f_damp_t = sqrt(one + rfac*(rfac - two) / bet_sqr);
-     orbit(2) = orbit(2) * f_damp_t;
-     orbit(4) = orbit(4) * f_damp_t;
-     orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
+      if(elrad .eq. zero) then 
+      call fort_warn('TWCPGO: ','Radiation effects ignored for solenoid '// &
+                         'with l=0, lrad=0 and radiate=true')
+    else
+       kx = ((sk**2)*orbit(1)-sk*orbit(4))*elrad;
+       ky = ((sk**2)*orbit(3)+sk*orbit(2))*elrad;
+       rfac = (arad * gamma**3 / three) * (kx**2 + ky**2) / elrad;
+       pt = orbit(6);
+       bet_sqr = (pt*pt + two*pt/bet0 + one) / (one/bet0 + pt)**2;
+       f_damp_t = sqrt(one + rfac*(rfac - two) / bet_sqr);
+       orbit(2) = orbit(2) * f_damp_t;
+       orbit(4) = orbit(4) * f_damp_t;
+       orbit(6) = orbit(6) * (one - rfac) - rfac / bet0;
+    endif
   endif
 
 end SUBROUTINE tmsol_th