BLCMDhelicaldipole.cc File Reference

#include <math.h>
#include "G4VisAttributes.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4VPhysicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4Color.hh"
#include "G4UserLimits.hh"
#include "G4Polymarker.hh"
#include "G4VVisManager.hh"
#include "BLElement.hh"
#include "BLElementField.hh"
#include "BLGlobalField.hh"
#include "BLParam.hh"
#include "BLManager.hh"
#include "BLHelicalUtils.hh"

Classes

class  BLCMDhelicaldipole
 BLCMDhelicaldipole implements helical dipole magnet with a cylindrical field volume. More...
class  HelicalDipoleField
 HelicalDipoleField represents one placement of a helicaldipole magnet. More...

Defines

#define _USE_MATH_DEFINES
#define MARKER_SIZE   5
#define cylPHI   0
#define cylRHO   1
#define cylZ   2
#define rectX   0
#define rectY   1
#define rectZ   2
#define ESSENTIALLY_ZERO   1.E-33
#define LEFT_HANDED_THREAD   (-1)
#define RIGHT_HANDED_THREAD   (+1)
#define MuonsInc_PITCH_CONVENTION   RIGHT_HANDED_THREAD
#define MuonsInc_Target_Radius_mm   159.15494
#define KBB_ENABLE   1
#define KBB_DISABLE   0
#define KBB_bugz   KBB_DISABLE

Functions

G4double modified_bessel0 (G4double x)
G4double modified_bessel1 (G4double x)
G4double modified_besselN (int o, G4double x)
G4double modified_besselN_x (int n, G4double x)
G4double derivative_modified_besselN (int n, G4double x)
G4ThreeVector DIPOLEFIELD (G4double bd, G4double k, G4double rho, G4double psiangle)
G4ThreeVector QUADRUPOLEFIELD (G4double bprime, G4double k, G4double rho, G4double psiangle)

Variables

const G4double PI = 3.141592653589793238462643
BLCMDhelicaldipole defaultHelicalDipole

Define Documentation

#define _USE_MATH_DEFINES

#define MARKER_SIZE   5

Referenced by BLCMDhelicalharmonic::EndOfRunAction(), and BLCMDhelicaldipole::EndOfRunAction().

#define cylPHI   0

#define cylRHO   1

#define cylZ   2

#define rectX   0

#define rectY   1

#define rectZ   2

#define ESSENTIALLY_ZERO   1.E-33

#define LEFT_HANDED_THREAD   (-1)

#define RIGHT_HANDED_THREAD   (+1)

#define MuonsInc_PITCH_CONVENTION   RIGHT_HANDED_THREAD

#define MuonsInc_Target_Radius_mm   159.15494

#define KBB_ENABLE   1

#define KBB_DISABLE   0

#define KBB_bugz   KBB_DISABLE

Referenced by HelicalDipoleField::addFieldValue(), and QUADRUPOLEFIELD().

Function Documentation

G4double modified_bessel0 ( G4double  x  ) 

Referenced by modified_besselN().

00407 { 
00408   /* Returns the modified Bessel function I_0(x) for any real x */
00409   /* modified bessel function of order 0 */ 
00410   /* slightly modified DBESI0 for type double, KBB 4/04  */
00411     double ax;
00412     double anser;
00413     double y;
00414  
00415     if ((ax=fabs(x)) < 3.75) {
00416       y = x/3.75;
00417       y*=y;
00418       anser = 1+y*(3.5156229+y*(3.0899424+y*(1.2067492
00419                +y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
00420     }
00421     else {
00422       y=3.75/ax;
00423       anser = (exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
00424               +y*(0.225319e-2+y*(-.157565e-2+y*(0.916281e-2
00425               +y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
00426               +y*0.392377e-2)))))))); 
00427     }
00428     return anser;
00429 }

G4double modified_bessel1 ( G4double  x  ) 

Referenced by modified_besselN().

00433 {   
00434   /* Returns the modified Bessel function I_1(x) for any real x */
00435   /* modified bessel function of order 1 */ 
00436   /* slightly modified DBESI1 for type double, KBB 4/04 */
00437     double ax;
00438     double anser;
00439     double y;
00440  
00441     if ((ax=fabs(x)) < 3.75) {
00442       y = x/3.75;
00443       y*=y;
00444       anser = ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
00445               +y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
00446     }
00447     else {
00448       y=3.75/ax;
00449       anser = 0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
00450               -y*0.420059e-2));
00451       anser = 0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
00452               +y*(0.163801e-2+y*(-0.1031555e-1+y*anser))));
00453       anser *= (exp(ax)/sqrt(ax));
00454          }
00455     return x < 0.0 ? -anser : anser;
00456 }

G4double modified_besselN ( int  o,
G4double  x 
)

References modified_bessel0(), modified_bessel1(), and modified_besselN_x().

Referenced by derivative_modified_besselN(), DIPOLEFIELD(), modified_besselN_x(), and QUADRUPOLEFIELD().

00459 {
00460   const char *CoDE="modified_besselN";
00461   int n;
00462   G4double modified_besselN_x( int n, G4double x ); /*prototype*/
00463   G4double Io=0;
00464 
00465   n= o-1;
00466   if( o>1 ) 
00467       Io= -2*n * modified_besselN_x(n,x) + modified_besselN(n-1,x);
00468   else if( o==1 )
00469       Io= modified_bessel1(x);
00470   else if( o==0 )
00471       Io= modified_bessel0(x);
00472   else
00473       printf("%s: bad order=%i\n", CoDE, n);
00474 
00475   return( Io );
00476 }

G4double modified_besselN_x ( int  n,
G4double  x 
)

References ESSENTIALLY_ZERO, and modified_besselN().

Referenced by derivative_modified_besselN(), DIPOLEFIELD(), modified_besselN(), and QUADRUPOLEFIELD().

00479 {
00480   G4double In_t;
00481   int i;
00482 
00483   if( -ESSENTIALLY_ZERO<x && x<ESSENTIALLY_ZERO )
00484     { 
00485       for(In_t=1/2.,i=2; i<=n; ++i)  In_t*= (x/2)/i;
00486     }
00487   else
00488     {
00489       In_t= modified_besselN(n,x)/x;
00490     }
00491 
00492   return( In_t );
00493 }

G4double derivative_modified_besselN ( int  n,
G4double  x 
)

References modified_besselN(), and modified_besselN_x().

00497 {
00498   const char *CoDE="derivative_modified_besselN";
00499 
00500   if( n<1 )  printf("%s: order#n < 1\n", CoDE);
00501   return( modified_besselN(n-1,x) - n * modified_besselN_x(n,x) );
00502 }

G4ThreeVector DIPOLEFIELD ( G4double  bd,
G4double  k,
G4double  rho,
G4double  psiangle 
)

References cylPHI, cylRHO, cylZ, modified_besselN(), and modified_besselN_x().

Referenced by HelicalDipoleField::addFieldValue().

00508 {
00509           const int PHI=cylPHI,RHO=cylRHO,Z=cylZ;
00510           G4double t,I0,I1_t,psi;
00511           G4ThreeVector Bcyl;
00512 
00513           t= k * rho;
00514 
00515           I0= modified_besselN(0,t);
00516 
00517           I1_t= modified_besselN_x(1,t);
00518 
00519           /* Helical Cooling Channel Simnulation Parameters and Fields */
00520           /* 6D Cooling Note                                           */
00521           /* using cylindrical components  - equation 1.2              */ 
00522 
00523           Bcyl[PHI]= 2 * bd * I1_t    * cos(psiangle);
00524 
00525           Bcyl[RHO]= 2 * bd * (I0-I1_t) * sin(psiangle);
00526 
00527           /* -2 * bd * I1(t) * cos(phiangle) <==> -t * Bcyl[PHI] */
00528 
00529           Bcyl[Z]=   - k * rho * Bcyl[PHI];
00530           //      printf("1; (%f, %f) %e %e %e \n",rho,psiangle,Bcyl[0],Bcyl[1],Bcyl[2]);
00531           return Bcyl;
00532 }

G4ThreeVector QUADRUPOLEFIELD ( G4double  bprime,
G4double  k,
G4double  rho,
G4double  psiangle 
)

References cylPHI, cylRHO, cylZ, KBB_bugz, modified_besselN(), modified_besselN_x(), and MuonsInc_Target_Radius_mm.

Referenced by HelicalDipoleField::addFieldValue().

00539 {
00540   /* 
00541    #note that the NUMERICAL value of bprime is NOT in [T/m];
00542    one must convert bprime/(tesla/meter) to get units of [T/m],
00543    but [T=0.001] and [meter=1000]! 
00544 
00545    corrected a mistake in the formula 4/27/04 KBB,RJ
00546   */
00547 
00548   const char *CoDE="QUADRUPOLEFIELD";
00549   const int PHI=cylPHI,RHO=cylRHO,Z=cylZ;
00550   static int messageCount=0;
00551   G4double t,I1,I2_t,Bprime;
00552   G4ThreeVector Bcyl;
00553 
00554 
00555   if( 0==messageCount++ )  printf("%s: version of 4/27/04\n", CoDE);
00556 
00557 
00558   //      <----- T/m--------->   <----- 1/1000 --->   <---- T/mm ------>
00559 
00560   Bprime= bprime/(tesla/meter) * (millimeter/meter) * (tesla/millimeter);
00561 
00562   t= 2 * k * rho;
00563   
00564   I2_t= modified_besselN_x(2,t);
00565     
00566   I1= modified_besselN(1,t);
00567   
00568   /* Helical Cooling Channel Simnulation Parameters and Fields */
00569   /* 6D Cooling Note                                           */
00570   /* using cylindrical components  - equation 1.7              */ 
00571   /* BUT THERE IS A TYPO in equantion 1.7 -  it should say     */
00572   /* 1/k NOT 1/2 in the 2nd line!                              */
00573   
00574   /* I (2*k*rho)/(k*k*rho)=> I (2*k*rho)/(2*k*rho) * 2/k */
00575   /*  2                       2                          */
00576  
00577   Bcyl[PHI] = 2/k * Bprime * I2_t * cos(2*psiangle); 
00578   
00579   Bcyl[RHO] = 1/k * Bprime * (I1 - 2*I2_t) * sin(2*psiangle);  // <== fixed! 
00580   
00581   Bcyl[Z] = - k * rho * Bcyl[PHI];
00582 
00583   if( KBB_bugz && 0.95*MuonsInc_Target_Radius_mm<rho && rho<1.05*MuonsInc_Target_Radius_mm ) 
00584     {
00585       printf("%s: bprime=%f[T/m] %f[native] t=%f I2_t=%f psiangle=%f ", 
00586              CoDE, bprime/(tesla/meter), bprime, t, I2_t, psiangle );
00587       printf("k=%f/mm rho=%fmm Bcyl=(%f,%f,%f)[native]=(%f,%f,%f)[T]\n", 
00588              k*mm, rho/mm, Bcyl[PHI], Bcyl[RHO], Bcyl[Z],
00589              Bcyl[PHI]/tesla, Bcyl[RHO]/tesla, Bcyl[Z]/tesla);
00590     }
00591 
00592   //  printf("2; (%f, %f) %e %e %e \n",rho,psiangle,Bcyl[0],Bcyl[1],Bcyl[2]);
00593   return Bcyl;
00594 }

Variable Documentation

const G4double PI = 3.141592653589793238462643

Referenced by HelicalDipoleField::addFieldValue().

BLCMDhelicaldipole defaultHelicalDipole

g4beamline