| |
Fundamental TechnologiesVoyager LECP Pages |
by Sheela Shodhan
*******************************************************************************
SUBROUTINE FDMOD(X,Y,Z,BX,BY,BZ)
* PURPOSE : THIS ROUTINE CALCULATES THE COMPONENTS OF THE MAGNETIC * * FIELD, Bx,By,Bz AT THE GIVEN POINT (X,Y,Z) IN SPACE. * * THIS FIELD IS PRODUCED BY TWO TILTED MAGNETS- ONE ROTATED * * CLOCKWISE, THE OTHER ROTATED COUNTER-CLOCKWISE BY THE SAME * * AMOUNT. * * A FIELD DUE TO A MAGNET LYING IN THE XY PLANE IS GIVEN BY: * * * * B(r) = 1 1 M(r') . n (r - r') da' * * --- 1 --------3 * * c 1 |r - r'| * * integral * * WHERE M(r') IS THE MAGNETISATION NORMAL TO THE SURFACE OF MAGNET. * * IN THIS ROUTINE, THE MAGNETISATION OF THE MAGNET IS LINEARLY * * FALLING. SO IT HAS THE FORM: M(r') = A(x'-x ) + B(y'- y ) + M * * o o o * * THEN, SUBSTITUTING THIS FORM IN THE ABOVE EXPRESSION YIELDS THE * * COMPONENTS OF THE MAGNETIC FIELD AT THE POINT (x,y,z). * * HOWEVER, SINCE THE FIELD HERE IS PRODUCED BY THE TWO TILTED MAGNETS, * * THIS POINT (x,y,z) IN THE EXPERIMENTAL SYSTEM HAS TO BE EXPRESSED * * IN THE COORDINATE SYSTEM THAT IS CENTRED AT EACH OF THE MAGNETS. * * THEREFORE, THIS POINT IS TRANSLATED AND ROTATED APPROPRIATELY INTO * * THE COORDINATE SYSTEMS OF EACH OF THE TWO MAGNETS. * * AFTER, THE FIELD DUE TO EACH OF THE MAGNETS IS COMPUTED, THE TOTAL * * FIELD IS DETERMINED BY ADDING IT UP DUE TO THE PRINCIPLE OF * * SUPERPOSITION. * * THE RELEVANT PARAMETERS:AMOUNT OF TRANSLATION,ROTATION, AND THE * * VALUES OF THE MAGNETISATION CONSTANTS A,B,AND M ARE INCLUDED IN * * o * * IN THE FILE FDMOD1.CMN. * * * * * * VARIABLES: * * INPUT: * * X,Y,Z : COORDINATES OF A POINT WHERE THE FIELD IS REQUIRED. * * OUTPUT: * * BX,BY,BZ : TOTAL FIELD COMPONENTS AT A POINT (X,Y,Z). * * OTHERS: * * ANG : ANGLE IN RADIANS BY WHICH THE MAGNETS ARE ROTATED. * * BX1,BX2,BXX1,BXX2,BXX3,BXX4,BXX5 : TEMPORARIES WHICH ARE USED TO STORE * * THE X-COMPONENT OF THE MAGNETIC FIELD. * * BY1,BY2,BYY1,BYY2,BYY3,BYY4,BYY5 : TEMPORARIES WHICH ARE USED TO STORE * * THE Y-COMPONENT OF THE MAGNETIC FIELD. * * BZ1,BZ2,BZZ1,BZZ2,BZZ3,BZZ4 : TEMPORARIES WHICH ARE USED TO STORE * * THE Z-COMPONENT OF THE MAGNETIC FIELD. * * XL,YL,ZL : ARRAYS WHOSE ELEMENTS CONTAIN THE AMOUNT BY WHICH THE * * CORRESPONDING PIECE OF THE MAGNET IS TRANSLATED. * *******************************************************************************
IMPLICIT NONE
REAL*8 XO,YO,ZO,X,Y,Z,Y1,Z1,Y2,Z2,YP,YN,XP,XN
REAL*8 BXX,BYY,BZZ,BX1,BX2,BY1,BY2,BZ1,BZ2
REAL*8 BX,BY,BZ
REAL*8 BXX1,BXX2,BXX3,BXX4,BXX5,BYY1,BYY2,BYY3,BYY4,BYY5
REAL*8 BZZ1,BZZ2,BZZ3,BZZ4
REAL*8 SQS,BS,COXZ,COYZ,DEN,NUM,TEMP,X1,X2
INTEGER I,LX,LY
INCLUDE 'FDMOD1.CMN'
C INITIALIZATION OF THE FIELD
BX=0.0D0
BY=0.0D0
BZ=0.0D0
DO LX=1,NX1
DO LY=1,NY1
C THE TRANSLATION OF THE COORDINATE SYSTEM
XO=X-XL(LX)
YO=Y+YL(LY)
ZO=Z+ZL(LY)
C THE ROTATION OF THE COORDINATE SYSTEM FOR LEFT MAGNET
X1=XO
Y1=YO*DCOS(ANG)-ZO*DSIN(ANG)
Z1=YO*DSIN(ANG)+ZO*DCOS(ANG)
C THE ROTATION OF THE COORDINATE SYSTEM FOR RIGHT MAGNET
X2=XO
Y2=YO*DCOS(ANG)+(ZO-SEPTS(LY))*DSIN(ANG)
Z2=-YO*DSIN(ANG)+(ZO-SEPTS(LY))*DCOS(ANG)
YP=Y1+WIDTHS(LY)
YN=Y1-WIDTHS(LY)
XP=XO+LENGTHS(LX)
XN=XO-LENGTHS(LX)
BZ1=BZZ(YP,XP,Z1)-BZZ(YP,XN,Z1)-BZZ(YN,XP,Z1)+BZZ(YN,XN,Z1)
BY1=BYY(YP,XN,Z1)-BYY(YP,XP,Z1)-BYY(YN,XN,Z1)+BYY(YN,XP,Z1)
BX1=BXX(YP,XN,Z1)-BXX(YP,XP,Z1)-BXX(YN,XN,Z1)+BXX(YN,XP,Z1)
d print *,'bx1 ',bx1,' by1 ',by1,' bz1 ',bz1
c to calculate the field due to linear
c correction in the magnetisation of the magnet.
BXX1 = B0 * ( (2.0d0*SQS(XN,Z1,Y1)) - (2.0d0*SQS(XP,Z1,Y1))
1 + SQS(XP,Z1,YN) + SQS(XP,Z1,YP)
1 - SQS(XN,Z1,YN) - SQS(XN,Z1,YP) )
BXX2 =( (2.0d0 * B0 * (Y1 - Y0)) *
1 (BXX(Y1,XP,Z1) - BXX(Y1,XN,Z1)) )
1 + (( B0 * (Y1 - Y0) ) *
1 ( BXX(YP,XN,Z1) + BXX(YN,XN,Z1)
1 - BXX(YP,XP,Z1) - BXX(YN,XP,Z1) ))
BXX3 = (( (A0/2.0d0) * (X0 - X1)) *
1 (BXX(YP,XP,Z1)
1 + BXX(YP,XN,Z1)
1 - BXX(YN,XP,Z1)
1 - BXX(YN,XN,Z1) ))
1 + (( A0*(X0-X1) ) *
1 ( BXX(YN,X1,Z1)
1 - BXX(YP,X1,Z1) ))
NUM = BS(YN,XP,Z1)*BS(YP,X1,Z1)*BS(YN,XN,Z1)
1 * BS(YP,X1,Z1)
DEN = BS(YN,X1,Z1)*BS(YP,XP,Z1)*BS(YN,X1,Z1)
1 * BS(YP,XN,Z1)
TEMP = NUM/DEN
BXX3 = BXX3 + ( ((A0*(X0-X1))/2.0) *dlog(TEMP) )
BXX4 = ( (A0 * (Y1 - WIDTHS(LY))) *
1 (BXX(XP,YN,Z1) + BXX(XN,YN,Z1)
1 -BXX(X1,YN,Z1) - BXX(X1,YN,Z1)) )
1 + ( (A0 * (Y1 + WIDTHS(LY))) *
1 (BXX(X1,YP,Z1) + BXX(X1,YP,Z1)
1 -BXX(XP,YP,Z1) - BXX(XN,YP,Z1)) )
COXZ = BZZ(X1,YN,Z1) + BZZ(X1,YN,Z1) + BZZ(XP,YP,Z1)
1 + BZZ(XN,YP,Z1) - BZZ(XN,YN,Z1)
1 - BZZ(XP,YN,Z1) - BZZ(X1,YP,Z1)
1 - BZZ(X1,YP,Z1)
BXX5 = A0 * Z1 * COXZ
d print *,'due to 1 magnet' d print *,'BXX1',BXX1,'BXX2',BXX2,'BXX3',BXX3 d print *,'BXX4',BXX4,'BXX5',BXX5 d print *,' '
c to compute the total x- field
BX1 = BX1 + BXX1 + BXX2 + BXX3 + BXX4 + BXX5
d print *,'due to 1 magnet' d print *,'BX1',BX1
c to compute the y-component of the field
BYY1 = A0 * ( (2.0d0*SQS(YN,Z1,X1)) - (2.0d0*SQS(YP,Z1,X1))
1 + SQS(YP,Z1,XN) + SQS(YP,Z1,XP)
1 - SQS(YN,Z1,XN) - SQS(YN,Z1,XP) )
BYY2 =( (2.0d0 * A0 * (X1 - X0)) *
1 (BXX(X1,YP,Z1) - BXX(X1,YN,Z1)) )
1 + (( A0 * (X1 - X0) ) *
1 ( BXX(XP,YN,Z1) + BXX(XN,YN,Z1)
1 - BXX(XP,YP,Z1) - BXX(XN,YP,Z1) ))
BYY3 = (( (B0/2.0d0) * (Y0 - Y1)) *
1 ( BXX(XP,YP,Z1)
1 + BXX(XP,YN,Z1)
1 - BXX(XN,YP,Z1)
1 - BXX(XN,YN,Z1) ))
1 + (( B0*(Y0-Y1) ) *
1 ( BXX(XN,Y1,Z1)
1 - BXX(XP,Y1,Z1) ))
NUM = BS(XN,YP,Z1)*BS(XN,YN,Z1)*BS(XP,Y1,Z1)
1 * BS(XP,Y1,Z1)
DEN = BS(XP,YP,Z1)*BS(XP,YN,Z1)*BS(XN,Y1,Z1)
1 * BS(XN,Y1,Z1)
TEMP = NUM/DEN
BYY3 = BYY3 + ( ((B0*(Y0-Y1))/2.0d0)*dlog(TEMP) )
BYY4 = ( (B0 * (X1 - LENGTHS(LX))) *
1 (BXX(YP,XN,Z1) + BXX(YN,XN,Z1)
1 -BXX(Y1,XN,Z1) - BXX(Y1,XN,Z1)) )
1 + ( (B0 * (X1 + LENGTHS(LX))) *
1 (BXX(Y1,XP,Z1) + BXX(Y1,XP,Z1)
1 -BXX(YP,XP,Z1) - BXX(YN,XP,Z1)) )
COYZ = BZZ(Y1,XN,Z1) + BZZ(Y1,XN,Z1) + BZZ(YP,XP,Z1)
1 + BZZ(YN,XP,Z1) - BZZ(YN,XN,Z1)
1 - BZZ(YP,XN,Z1) - BZZ(Y1,XP,Z1)
1 - BZZ(Y1,XP,Z1)
BYY5 = B0 * Z1 * COYZ
c to compute the total y-field
BY1 = BY1 + BYY1 + BYY2 + BYY3 + BYY4 + BYY5
d print *,'due to 1 magnet' d print *,'BY1',BY1,'BYY1',BYY1,'BYY2',BYY2 d print *,'BYY3',BYY3,'BYY4',BYY4,'BYY5',BYY5 d print *,' '
c next, to compute the z-component of the field
BZZ1 = ( (2.0d0*A0*Z1) *
1 (BXX(YN,X1,Z1) - BXX(YP,X1,Z1)) )
1 + ( (A0*Z1) *
1 (BXX(YP,XP,Z1) + BXX(YP,XN,Z1)
1 - BXX(YN,XP,Z1) - BXX(YN,XN,Z1)) )
BZZ2 = A0 * (X1 - X0) * COXZ
BZZ3 = ( (B0 * Z1) *
1 ( BXX(XN,Y1,Z1)
1 - BXX(XP,Y1,Z1) ) )
1 + ( ((B0 * Z1)/2.0d0) *
1 (BXX(XP,YP,Z1)
1 + BXX(XP,YN,Z1)
1 - BXX(XN,YP,Z1)
1 -BXX(XN,YN,Z1) ) )
NUM = BS(XP,Y1,Z1)*BS(XP,Y1,Z1)*BS(XN,YP,Z1)
1 * BS(XN,YN,Z1)
DEN = BS(XP,YP,Z1)*BS(XN,Y1,Z1)*BS(XN,Y1,Z1)
1 * BS(XP,YN,Z1)
TEMP = dlog(NUM/DEN)
TEMP = ((B0*Z1)/2.0d0) * TEMP
BZZ3 = BZZ3 + TEMP
BZZ4 = (B0* (Y0-Y1)) *
1 (BZZ(XN,YN,Z1) + BZZ(XP,Y1,Z1)
1 + BZZ(XP,Y1,Z1) + BZZ(XN,YP,Z1)
1 - BZZ(XN,Y1,Z1) - BZZ(XN,Y1,Z1)
1 - BZZ(XP,YN,Z1) - BZZ(XP,YP,Z1))
c to compute the total z-field
BZ1 = BZ1 + BZZ1 + BZZ2 + BZZ3 + BZZ4
d print *,'due to 1 magnet' d print *,'BZ1',BZ1,'BZZ1',BZZ1,'BZZ2',BZZ2 d print *,'BZZ3',BZZ3,'BZZ4',BZZ4 d print *,' '
c for the 2 magnet
YP=Y2+WIDTHS(LY)
YN=Y2-WIDTHS(LY)
BZ2=BZZ(YP,XP,Z2)-BZZ(YP,XN,Z2)-BZZ(YN,XP,Z2)+BZZ(YN,XN,Z2)
BY2=BYY(YP,XN,Z2)-BYY(YP,XP,Z2)-BYY(YN,XN,Z2)+BYY(YN,XP,Z2)
BX2=BXX(YP,XN,Z2)-BXX(YP,XP,Z2)-BXX(YN,XN,Z2)+BXX(YN,XP,Z2)
d print *,' bx2 ',bx2,' by2 ',by2,' bz2 ',bz2
c c add the correction due to the linear magnetization
BXX1 = B0 * ( (2.0d0*SQS(XN,Z2,Y2)) - (2.0d0*SQS(XP,Z2,Y2))
1 + SQS(XP,Z2,YN) + SQS(XP,Z2,YP)
1 - SQS(XN,Z2,YN) - SQS(XN,Z2,YP) )
BXX2 =( (2.0d0 * B0 * (Y2 - Y0)) *
1 (BXX(Y2,XP,Z2) - BXX(Y2,XN,Z2)) )
1 + (( B0 * (Y2 - Y0) ) *
1 ( BXX(YP,XN,Z2) + BXX(YN,XN,Z2)
1 - BXX(YP,XP,Z2) - BXX(YN,XP,Z2) ))
BXX3 = (( (A0/2.0d0) * (X0 - X2)) *
1 (BXX(YP,XP,Z2)
1 + BXX(YP,XN,Z2)
1 - BXX(YN,XP,Z2)
1 - BXX(YN,XN,Z2) ))
1 + (( A0*(X0-X2) ) *
1 ( BXX(YN,X2,Z2)
1 - BXX(YP,X2,Z2) ))
NUM = BS(YN,XP,Z2)*BS(YP,X2,Z2)*BS(YN,XN,Z2)
1 * BS(YP,X2,Z2)
DEN = BS(YN,X2,Z2)*BS(YP,XP,Z2)*BS(YN,X2,Z2)
1 * BS(YP,XN,Z2)
TEMP = NUM/DEN
BXX3 = BXX3 + ( ((A0*(X0-X2))/2.0d0) *dlog(TEMP) )
BXX4 = ( (A0 * (Y2 - WIDTHS(LY))) *
1 (BXX(XP,YN,Z2) + BXX(XN,YN,Z2)
1 -BXX(X2,YN,Z2) - BXX(X2,YN,Z2)) )
1 + ( (A0 * (Y2 + WIDTHS(LY))) *
1 (BXX(X2,YP,Z2) + BXX(X2,YP,Z2)
1 -BXX(XP,YP,Z2) - BXX(XN,YP,Z2)) )
COXZ = BZZ(X2,YN,Z2) + BZZ(X2,YN,Z2) + BZZ(XP,YP,Z2)
1 + BZZ(XN,YP,Z2) - BZZ(XN,YN,Z2)
1 - BZZ(XP,YN,Z2) - BZZ(X2,YP,Z2)
1 - BZZ(X2,YP,Z2)
BXX5 = A0 * Z2 * COXZ
c to compute the total x- field
BX2 = BX2 + BXX1 + BXX2 + BXX3 + BXX4 + BXX5
d print *,'due to 2 magnet' d print *,'BX2',BX2,'BXX1',BXX1,'BXX2',BXX2 d print *,'BXX3',BXX3,'BXX4',BXX4 d print *,' '
c to compute the y-component of the field
BYY1 = A0 * ( (2.0d0*SQS(YN,Z2,X2)) - (2.0d0*SQS(YP,Z2,X2))
1 + SQS(YP,Z2,XN) + SQS(YP,Z2,XP)
1 - SQS(YN,Z2,XN) - SQS(YN,Z2,XP) )
BYY2 =( (2.0d0 * A0 * (X2 - X0)) *
1 (BXX(X2,YP,Z2) - BXX(X2,YN,Z2)) )
1 + (( A0 * (X2 - X0) ) *
1 ( BXX(XP,YN,Z2) + BXX(XN,YN,Z2)
1 - BXX(XP,YP,Z2) - BXX(XN,YP,Z2) ))
BYY3 = (( (B0/2.0d0) * (Y0 - Y2)) *
1 ( BXX(XP,YP,Z2)
1 + BXX(XP,YN,Z2)
1 - BXX(XN,YP,Z2)
1 - BXX(XN,YN,Z2) ))
1 + (( B0*(Y0-Y2) ) *
1 ( BXX(XN,Y2,Z2)
1 - BXX(XP,Y2,Z2) ))
NUM = BS(XN,YP,Z2)*BS(XN,YN,Z2)*BS(XP,Y2,Z2)
1 * BS(XP,Y2,Z2)
DEN = BS(XP,YP,Z2)*BS(XP,YN,Z2)*BS(XN,Y2,Z2)
1 * BS(XN,Y2,Z2)
TEMP = NUM/DEN
BYY3 = BYY3 + ( ((B0*(Y0-Y2))/2.0)*dlog(TEMP) )
BYY4 = ( (B0 * (X2 - LENGTHS(LX))) *
1 (BXX(YP,XN,Z2) + BXX(YN,XN,Z2)
1 -BXX(Y2,XN,Z2) - BXX(Y2,XN,Z2)) )
1 + ( (B0 * (X2 + LENGTHS(LX))) *
1 (BXX(Y2,XP,Z2) + BXX(Y2,XP,Z2)
1 -BXX(YP,XP,Z2) - BXX(YN,XP,Z2)) )
COYZ = BZZ(Y2,XN,Z2) + BZZ(Y2,XN,Z2) + BZZ(YP,XP,Z2)
1 + BZZ(YN,XP,Z2) - BZZ(YN,XN,Z2)
1 - BZZ(YP,XN,Z2) - BZZ(Y2,XP,Z2)
1 - BZZ(Y2,XP,Z2)
BYY5 = B0 * Z2 * COYZ
c to compute the total y-field
BY2 = BY2 + BYY1 + BYY2 + BYY3 + BYY4 + BYY5
d print *,'due to 2 magnet' d print *,'BY2',BY2,'BYY1',BYY1,'BYY2',BYY2 d print *,'BYY3',BYY3,'BYY4',BYY4 d print *,' '
c next, to compute the z-component of the field
BZZ1 = ( (2.0d0*A0*Z2) *
1 (BXX(YN,X2,Z2) - BXX(YP,X2,Z2)) )
1 + ( (A0*Z2) *
1 (BXX(YP,XP,Z2) + BXX(YP,XN,Z2)
1 - BXX(YN,XP,Z2) - BXX(YN,XN,Z2)) )
BZZ2 = A0 * (X2 - X0) * COXZ
BZZ3 = ( (B0 * Z2) *
1 ( BXX(XN,Y2,Z2)
1 - BXX(XP,Y2,Z2) ) )
1 + ( ((B0 * Z2)/2.0d0) *
1 (BXX(XP,YP,Z2)
1 + BXX(XP,YN,Z2)
1 - BXX(XN,YP,Z2)
1 -BXX(XN,YN,Z2) ) )
NUM = BS(XP,Y2,Z2)*BS(XP,Y2,Z2)*BS(XN,YP,Z2)
1 * BS(XN,YN,Z2)
DEN = BS(XP,YP,Z2)*BS(XN,Y2,Z2)*BS(XN,Y2,Z2)
1 * BS(XP,YN,Z2)
TEMP = dlog(NUM/DEN)
TEMP = ((B0*Z2)/2.0d0) * TEMP
BZZ3 = BZZ3 + TEMP
BZZ4 = (B0* (Y0-Y2)) *
1 (BZZ(XN,YN,Z2) + BZZ(XP,Y2,Z2)
1 + BZZ(XP,Y2,Z2) + BZZ(XN,YP,Z2)
1 - BZZ(XN,Y2,Z2) - BZZ(XN,Y2,Z2)
1 - BZZ(XP,YN,Z2) - BZZ(XP,YP,Z2))
c to compute the total z-field
BZ2 = BZ2 + BZZ1 + BZZ2 + BZZ3 + BZZ4
d print *,'due to 2 magnet' d print *,'BZ2',BZ2,'BZZ1',BZZ1,'BZZ2',BZZ2 d print *,'BZZ3',BZZ3,'BZZ4',BZZ4 d print *,' '
BZ = BZ + (MAGS(LX,LY) * (BZ2-BZ1))
BY = BY + (MAGS(LX,LY) * (BY2-BY1))
BX = BX + (MAGS(LX,LY) * (BX2-BX1))
D PRINT *,'MAGS(',LX,LY,')',MAGS(LX,LY)
D PRINT *,' X ',X,' Y ',Y,' Z ',Z
D PRINT *,'BZ1:',BZ1,' BZ2 ',BZ2,' BZ ',BZ
D PRINT *,'BY1:',BY1,' BY2 ',BY2,' BY ',BY
D PRINT *,'BX1:',BX1,' BX2 ',BX2,' BX ',BX
END DO
END DO
RETURN
END
C---------------------------------------------------------------------
REAL*8 FUNCTION BZZ(P1,P2,P3)
REAL*8 TEMP,P1,P2,P3
IF (P3 .EQ. 0.0D0) THEN
TEMP = PI/2.0D0
ELSE
TEMP = (P1 * P2)/(P3 * DSQRT(P1**2+P2**2+P3**2))
TEMP = DATAN(TEMP)
END IF
BZZ=TEMP
RETURN
END
C---------------------------------------------------------------------
c the function SQS which computes the c square root of the sum of the squares of the input c numbers.
REAL*8 function SQS(P1,P2,P3)
REAL*8 P1,P2,P3
SQS = DSQRT( (P1*P1) + (P2*P2) + (P3*P3) )
return
end
C---------------------------------------------------------------------
c this is the function that calculates c (p1 - DSQRT(p1*p1 + p2*p2 + p3*p3))
REAL*8 function BS(P1,P2,P3)
REAL*8 P1,P2,P3,TEMP
TEMP = DSQRT((P1*P1) + (P2*P2) + (P3*P3))
BS = P1 - TEMP
return
end
C---------------------------------------------------------------------
REAL*8 FUNCTION BXX(P1,P2,P3)
REAL*8 TEMP,P1,P2,P3
TEMP=DSQRT(P1**2+P2**2+P3**2)
BXX=DLOG(TEMP+P1)
RETURN
END
C--------------------------------------------------------------------
REAL*8 FUNCTION BYY(P1,P2,P3)
REAL*8 TEMP,P1,P2,P3
TEMP=DSQRT(P1**2+P2**2+P3**2)
BYY=DLOG(TEMP+P2)
RETURN
END
C--------------------------------------------------------------------
C--------------------------------------------------------------------
Return to thesis table of contents.
Return to Voyager LECP Data Analysis Handbook Table of
Contents.
Return to Fundamental Technologies Home Page.
Last modified 12/9/02, Tizby Hunt-Ward
tizby@ftecs.com