Hi
When I try to compute the Gauss formula to calculate the number of turn's wraps between two solenoids. if I don't put spaces between the turns (here a linear function,t(turnSpace=0), then the Integer value is inside the error range of the monte_carlo_integral() function. But if I introduce a space between the turns (turnSpace=10242*pi), then appears an error in the calculation of the Monte Carlo integral value which is out of the error range !. (comment lines 61 and 62 to introduce a space between the turns of both coils) Where is the error ? in my computation (most likely !)? in the monte_carlo_integral function ?
with another software when integrating the function below between t=[0,12pi] and s=[0,14pi] do we find the integer value 42=6*7 +- Error with the integer value in the error range?
print("https://www.maths.ed.ac.uk/~v1ranick/papers/ricca.pdf")
allVars=var('alpha,beta,gamma,'+
't,t0,t1,nt,shiftXt,shiftYt,shiftZt,Rt,turnSpaceT,'+
's,s0,s1,ns,shiftXs,shiftYs,shiftZs,Rs,turnSpaceS',domain='real')
assume(nt>0) ;assume(ns >0);
precis=50
#########################################
def addAxis(plt,textToPlt,xBound,yBound,zBound) :
# add axis for real part
plt += arrow3d((0, 0, 0), (xBound, 0, 0), color='blue')
plt += text3d("X", (xBound, 0, 0), color='blue')
plt += arrow3d((0, 0, 0), (0, yBound, 0), color='green')
plt += text3d("Y", (0, yBound, 0), color='green')
plt += arrow3d((0, 0, 0), (0, 0, zBound), color='red')
plt += text3d("Z", (0,0, zBound), color='red')
plt+=textToPlt
return plt
def RotateVector(V,anglesDic) :
#rotation angle alpha around X axis
RotalphaX=matrix(SR,[[1,0,0],[0,cos(alpha),-sin(alpha)],[0,sin(alpha),cos(alpha)]]).subs(anglesDic)
#rotation angle beta around Y axis
RotbetaY=matrix(SR,[[cos(beta),0,sin(beta)],[0,1,0],[-sin(beta),0,cos(beta)]]).subs(anglesDic)
#rotation angle gamma around Z axis
RotgammaZ=matrix(SR,[[cos(gamma),-sin(gamma),0],[sin(gamma),cos(gamma),0],[0,0,1]]).subs(anglesDic)
Vx= (RotalphaX*V)
Vxy= (RotbetaY*Vx)
Vxyz= (RotgammaZ*Vxy)
return vector(Vxyz)
anglesDict={alpha:0,beta:0,gamma:0}
anglesDics={alpha:pi/4,beta:0,gamma:pi/6}
#########################################
xyBound=1.2 ;
thck=8
gammat0 = vector(SR,[shiftXt+Rt*cos(t),shiftYt+Rt*sin(t),t*turnSpaceT]) # in xy plane
gammat=RotateVector(gammat0,anglesDict)
gammas0 = vector(SR,[shiftXs+Rs*cos(s),shiftYs+Rs*sin(s),s*turnSpaceS])
gammas= RotateVector(gammas0,anglesDics)
show("gammat : \t " ,gammat)
show("gammas : \t " ,gammas)
# set here below the number n of coil turns of the blue solenoid)
# the radius Rt and the turns space turnSpaceT
tDic={nt:6 , Rt:1 , turnSpaceT:2*pi/512 , shiftXt:1/2 , shiftYt:0 , shiftZt:0}
# change here below the number n of coil turns of the black solenoid)
# the radius Rs and the turns space turnSpaceS
sDic={ns:7 , Rs:1 , turnSpaceS:2*pi/512 , shiftXs:0 , shiftYs:0 ,shiftZs:0}
# no space between spires then Monte Carlo integration OK
#tDic={nt:6 , Rt:1 , turnSpaceT:0 , shiftXt:1/2 , shiftYt:0 , shiftZt:0}
#sDic={ns:7 , Rs:1 , turnSpaceS:0 , shiftXs:0 , shiftYs:0 ,shiftZs:0}
show("gammat with coil's values : \t " ,gammat.subs(tDic))
show("gammas with coil's values : \t " ,gammas.subs(sDic))
# set the circles centers
centerGammat=RotateVector(vector([shiftXt,shiftYt,shiftZt]).subs(tDic),anglesDict)
centerGammas=RotateVector(vector([shiftXs,shiftYs,shiftZs]).subs(sDic),anglesDics)
textToPlt = text3d("Gauss turns number", (xyBound, xyBound, xyBound), fontsize='300%', fontweight=800)
plt=parametric_plot3d(gammat.subs(tDic), (0,nt.subs(tDic)*2*pi),color="blue",thickness=thck)
plt+=parametric_plot3d(gammas.subs(sDic), (0,ns.subs(sDic)*2*pi),color="grey",thickness=thck)
plt+=addAxis(plt,textToPlt,1,1,1)
plt+=plot(centerGammat,color="blue")
plt+=point3d(centerGammat,size=20,color="blue")
plt+=plot(centerGammas,color="black")
plt+=point3d(centerGammas,size=20,color="black")
plt.show(viewer='threejs')
integrationFactorFunction=(1/(4*pi*((gammas-gammat).norm())^3)).subs(tDic).subs(sDic)
# here choose the integration function factor (integrationFactorFunction0 or integrationFactorFunction1)
M=matrix([ [ gammat.diff(t)[0],gammat.diff(t)[1],gammat.diff(t)[2] ],
[ gammas.diff(s)[0],gammas.diff(s)[1],gammas.diff(s)[2] ],
[(gammat[0]-gammas[0]), (gammat[1]-gammas[1]),(gammat[2]-gammas[2]) ] ])
precis=40 # precision approximation
# solenoid with space turn length
show(LatexExpr(r"tLen=\int^{b}_{a} \sqrt{ x_t′(t)^2 + y_t′(t)^2 + z_t′(t)^2 } dt \, \, \," +
r"sLen=\int^{b}_{a} \sqrt{ x_s′(t)^2 + y_s′(t)^2 + z_s′(t)^2 } dt "))
show(" ################ compute solenoid length ########################")
tLenDt=sqrt((M[0][0]^2 + M[0][1]^2 + M[0][2]^2).subs(tDic))
show(tLenDt)
solenoidTLen=integrate(tLenDt, (t, 0, nt.subs(tDic)*2*pi ) )
show("solenoid T length : \t ",solenoidTLen, " \t numerical value : \t ",solenoidTLen.n(precis))
sLenDt=sqrt((M[1][0]^2 + M[1][1]^2 + M[1][2]^2).subs(sDic))
show(sLenDt)
solenoidSLen=integrate(sLenDt, (s, 0, ns.subs(sDic)*2*pi ) )
show("solenoid S length : \t ",solenoidSLen, " \t numerical value : \t ",solenoidSLen.n(precis))
show("matrix M :")
show(M)
turnsNumber=integrationFactorFunction * M.det()
show("##################### turnsNumber #####################")
turnsNumberNumDerivative=(turnsNumber.subs(tDic)).subs(sDic)
textToPlts = text3d("surface to sum", (xyBound, xyBound, xyBound), fontsize='300%', fontweight=800,transparency=0.5)
var('x,y')
show("turns number fonction to integrate :",turnsNumberNumDerivative)
show(turnsNumberNumDerivative)
show("rough monte_carlo integral approximation : \t",
monte_carlo_integral(turnsNumberNumDerivative(x=t,y=s), [0,0], [ nt.subs(tDic)*2*pi, ns.subs(sDic)*2*pi], 10^8))
#for algo in ['plain', 'miser', 'vegas']: # abs tol 0.01
# show("rough monte_carlo integral approximation : \t",algo,
# monte_carlo_integral(turnsNumberNumDerivative(x=t,y=s), [0,0], [ nt.subs(tDic)*2*pi, ns.subs(sDic)*2*pi],
# 10^6, algorithm=algo))
show(r"""When I try to compute the Gaussian formula to calculate the number of turn's wraps between two solenoids.
if I don't put spaces between the turns (here a linear function,t*(turnSpace=0),
then the Integer value is inside the error range of the monte_carlo_integral() function.
But when I introduce a space between the turns(turnSpace=1024*2*pi), then appears an error in the computation
of the Monte Carlo integral value which is out of the error range !.
(comment lines 61 and 62 to introduce a space between the turns of both coils)
Where is the error ? in my computation ? in the monte_carlo_integral function ? in the Gauss formula ?""")
