{----------------------------------------------------------- Program to calculate the B-field from a solenoid (1 layer). The solenoid is placed in a (x,y,z) coordinate system with the center axis along the x-axis and with the center at origo. Current is flowing in positive direction of the x-axis and turns are wound clockwise in this direction. ------------------------------------------------------------} Program Solenoid; uses crt; { Edit constants before compilation/execution or rewrite the below to a question/answer dialog } const length: extended = 1; { Length of coil [m] } turns: longint = 1000; { Number of turns } diam: extended = 0.01; { coil diameter [m] } Bx: extended = 0; { (x,y,z) position in where the } By: extended = 0; { B-field is calculated [m] } Bz: extended = 0; I: extended = 250; { Current [A] } type vect3 = array[1..3] of extended; { (x,y,z) array } var Teta,dTeta: extended; posx: extended; { Position on x-axis } dx: extended; { stepsize on x-axis } B: vect3; { calculated B-field vector } Bpos: vect3; { Position of calculated B-field vector } dI: vect3; { I*ds vector } Ipos: vect3; { Position of dI vector } steps: longint; { Number of steps during integration } j: longint; { loop counter } { Data initialization } Procedure Init; begin steps := turns*360; { stepsize in windings = 1 deg. } dx := length/steps; dTeta := 2*pi/360; Teta := 0; { Actual angle for integration in windings } posx := -0.5*length; { x-position for negative end of coil } Bpos[1] := Bx; Bpos[2] := By; Bpos[3] := Bz; B[1] := 0; B[2] := 0; B[3] := 0; end; { Cross product } Procedure xprod( var a,b,c: vect3 ); { c = a x b } begin c[1] := a[2]*b[3] - a[3]*b[2]; c[2] := a[3]*b[1] - a[1]*b[3]; c[3] := a[1]*b[2] - a[2]*b[1]; end; { Calculate 3D distance } Function dist3D( var a,b: vect3 ): extended; begin dist3D := sqrt( sqr( abs( a[1]-b[1] )) + sqr( abs( a[2]-b[2] )) + sqr( abs( a[3]-b[3] )) ); end; { Calculate dI = I * ds } Procedure dI3D( I,alpha: extended; var dI: vect3 ); begin dI[1] := dx*I; dI[2] := -0.5*diam*sin(alpha)*dTeta*I; dI[3] := 0.5*diam*cos(alpha)*dTeta*I; end; { Increment x-position and Teta-angle } Procedure StepInteg; begin Teta := Teta + dTeta; if (Teta>2*pi) then Teta := Teta - 2*pi; posx := posx + dx; end; { Do integration step } Procedure IntegB; var r,r3: extended; rvect,dB: vect3; begin dI3D( I,Teta,dI ); { Construct dI-vector } Ipos[1] := posx; Ipos[2] := 0.5*diam*cos(Teta); { Calculate position of I-vector } Ipos[3] := 0.5*diam*sin(Teta); rvect[1] := Bpos[1] - Ipos[1]; rvect[2] := Bpos[2] - Ipos[2]; { Calculate r-vector } rvect[3] := Bpos[3] - Ipos[3]; r := dist3D(Bpos,Ipos); { Calculate | r | } r3 := r*r*r; xprod(dI,rvect,dB); { dB = Cross product } B[1] := B[1] + dB[1]*1E-7/r3; B[2] := B[2] + dB[2]*1E-7/r3; { Accumulate B-vector } B[3] := B[3] + dB[3]*1E-7/r3; end; { MAIN } begin Init; clrscr; gotoxy(1,10); write('Calculating...... '); for j := 1 to steps do begin IntegB; { Integration loop } StepInteg; end; gotoxy(1,10); writeln('pos-x[m]':12,'Bx[H]':15,'By[H]':15,'Bz[H]':15); writeln; writeln(posx:15:10,B[1]:15:10,B[2]:15:10,B[3]:15:10); writeln; writeln('':5,'Done'); readln; { program ends when 'enter' key pressed } end.