ULYSSES
Ulysses HISCALE Data Analysis Handbook
Appendix 9 Geometric Factor Study for the Deflected and Unscattered Electrons of HISCALE (Buckley MS Thesis)
A9.13 Appendix: Computer Programs
* PROGRAM I.3 *
******************************************************************************* * FDMOD1 * * THE MAIN PURPOSE OF THIS SUBROUTINE IS TO COMPUTE THE MAGNETIC FIELD AT * * A PARTICULAR LOCATION -- THIS LOCATION IS READ IN AS A COORDINATE. * * AGAIN, SEE SHODHAN'S THESIS FOR FURTHER COMMENTS. * *******************************************************************************
SUBROUTINE FDMOD(X,Y,Z,BX,BY,BZ)
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 '[HONG.ulybfield]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 D print *,A0,B0,MAGS(1,1)
c*********************************************************************
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*********************************************************************
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
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 C
IF (P3 .EQ. 0.0D0) THEN temp = p1/2.0D0 c 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 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************************************************************** C REAL*8 FUNCTION BXX(P1,P2,P3)
REAL*8 TEMP,P1,P2,P3 C TEMP=DSQRT(P1**2+P2**2+P3**2)
BXX=DLOG(TEMP+P1) RETURN END C C ****************************************************************** C REAL*8 FUNCTION BYY(P1,P2,P3) REAL*8 TEMP,P1,P2,P3 C TEMP=DSQRT(P1**2+P2**2+P3**2)
BYY=DLOG(TEMP+P2) RETURN END C ****************************************************************
Return to the Table of Contents for Buckley's Thesis
Return to HISCALE List of Appendices
Return to Ulysses HISCALE Data Analysis Handbook Table of Contents
Updated 8/8/19, Cameron Crane
QUICK FACTS
Manufacturer:
ESA provided the Ulysses spacecraft, NASA provided the power
supply, and various others provided its instruments.
Mission End Date: June 30, 2009
Destination: The inner heliosphere of the sun away from the ecliptic plane
Orbit: Elliptical orbit transversing the polar regions of the sun outside of the ecliptic plane
Mission End Date: June 30, 2009
Destination: The inner heliosphere of the sun away from the ecliptic plane
Orbit: Elliptical orbit transversing the polar regions of the sun outside of the ecliptic plane