Okay. Here you go. MonteHelix is the main routine which calls other routines (also shown).
There are some calls to some custom random number generator routines since I’ve found that the built-in random functions are not up to the task. If you want to see those, I can post them too.
[code]Function MonteHelix(n As int64,pitch As Double,r As Double,dw As Double,Maxerr As Double, Nmonte As int64) As double
'Monte carlo variables
dim sx,sy,sye,tx,ty,tye,montesumC,montesumU,ra,rb,y,yc,psi,cospsiinv,rra,rrb,p,LeftTermC,LeftTermU,L,Lc As Double
dim cIterTerm, uIterTerm, cSqSum,uSqSum As Double
dim ms As Double
dim iMonte,nupdate As Int64
dim tab As String
tab=chr(9)
nupdate=max(min(10000,Nmonte\1000),10)
SetRndSeed 'initialize random number generator
'calc pitch angle adjustments
psi= ATan(pitch/(pi2r))
cospsiinv=1/cos(psi)
OutText=OutText+EndOfLine+“psi factor: “+Format(cospsiinv,“0.0000000000”)
'outer montecarlo random filament generator loop
montesumU=0
montesumC=0
cSqSum=0
uSqSum=0
ms=Microseconds
for iMonte=1 to Nmonte
'Generate two random test points S and T
do
tx=rnd2-0.5
ty=rnd2-0.5
tye=sqrt(0.25-tx^2)
loop until (ty<=tye and ty>=-tye)
do
sx=rnd2-0.5
sy=rnd2-0.5
sye=sqrt(0.25-sx^2)
Loop Until (sy<=sye and sy>=-sye)
'calculate filament offset
y=(tx-sx)dw
yc=ycospsiinv 'psi corrected version of axial offset
ra=r+sydw
rb=r+tydw
L=Msnow(n,pitch,ra,rb,y,MaxErr)
'Lc=Msnow(n,pitch,ra,rb,yc,MaxErr)
rra=rara
rrb=rbrb
p=pitch/2/pi
LeftTermC=1e-7sqrt(rra+pp)sqrt(rrb+pp)/sqrt((1+pp/rra)(1+pp/rrb))
LeftTermU=1e-7sqrt(rrarrb)/sqrt((1+pp/rra)(1+pp/rrb))
cIterTerm=LeftTermCL
uIterTerm=LeftTermUL
montesumC=montesumC+cIterTerm
montesumU=montesumU+uIterTerm
cSqSum=cSqSum+cIterTermcIterTerm
uSqSum=uSqSum+uIterTermuIterTerm
if iMonte mod nupdate=0 then
OutText=OutText+EndOfLine+str(iMonte)+”, “+Format(montesumC/iMonte,”-#.############e”)+“, “+Format(montesumU/iMonte,”-#.############e”)
OutText=OutText+“, “+Format(cSqSum/iMonte,”-#.############e”)+“, “+Format(uSqSum/iMonte,”-#.############e”)
OutText=OutText+ “, “+ElapsedTime(ms)
end if
next
OutText=OutText+ EndOfLine+str(Nmonte)+tab+Format(montesumC/Nmonte,”-#.############e”)+tab+Format(montesumU/Nmonte,“-#.############e”)
OutText=OutText+ tab+Format(cSqSum/Nmonte,“-#.############e”)+tab+Format(uSqSum/Nmonte,“-#.############e”)
return montesumC/Nmonte
End Function
Function Msnow(n As Double, pitch As Double, ra As Double, rb As Double, y As Double, Maxerr As Double) As Double
’ Uses xxxxx original xxxxx formula to calculate xxxxxx between
’ two xxxxxx of arbitrary radii and axial offset
’ with xxxxx offset radially and axially according to optimum xxxxx rotation angle
’ evaluated using Simpson’s rule, and xxxxxx
’ This is the inner function of the full integration
Dim nstr as string
Dim p,Ntp,rra,rrb,rab,a,b,grandtotal,dx As Double
Dim CurrentErr,Sum2,LastIntg,Sum,Integral,LeftTermC,LeftTermU As Double
Dim m,mtot,i As Int64
p=pitch/2/pi
Ntp=2piN
nstr=“”
rra=rara
rrb=rbrb
rab=rarb
LeftTermC=1e-7sqrt(rra+pp)sqrt(rrb+pp)/sqrt((1+pp/rra)(1+pp/rrb))
LeftTermU=1e-7sqrt(rrarrb)/sqrt((1+pp/rra)(1+pp/rrb))
//Limits of integration are a (lower) and b (upper)
//The region 0…pi/2 is evaluated first because it may need the most subdivisions
a=0
b=pi/2
If b>Ntp then b=Ntp //ensure max value for b is Ntp
grandtotal=0
mtot=0
while a<Ntp // ensure ‘a’ never goes past upper limit
dx=b-a //range of integration
m=1
CurrentErr=2MaxErr
Sum2=HlxIntgrndRad(1,-a,Ntp,p,rra,rrb,rab,y)+HlxIntgrndRad(1,-b,Ntp,p,rra,rrb,rab,y)+HlxIntgrndRad(-1,a,Ntp,p,rra,rrb,rab,y)+HlxIntgrndRad(-1,b,Ntp,p,rra,rrb,rab,y)
'Initialize LastResult to trapezoidal area for test purposes
LastIntg=Sum2/2dx
While CurrentErr>MaxErr or m<512
m=2m
dx=dx/2
Sum=0
for i=1 to m step 2
Sum=Sum+HlxIntgrndRad(1,-idx-a,Ntp,p,rra,rrb,rab,y)+HlxIntgrndRad(-1,idx+a,Ntp,p,rra,rrb,rab,y)
next
Integral=(4Sum+Sum2)dx/3
CurrentErr=Abs((Integral)/(LastIntg)-1)
LastIntg=Integral
Sum2=Sum2+Sum2
wend
grandtotal=grandtotal+Integral
a=b
b=b2
If b>Ntp then b=Ntp
'nstr=nstr+str(m)+" "
mtot=mtot+m+1
wend
return grandtotal
End Function
Function HlxIntgrndRad(sn As Integer,phi As Double,N as Double,p As Double,rra As Double,rrb As Double,rab As Double,y As Double) As Double
return (N+snphi)(cos(phi)+pp/rab)/sqrt(rra+rrb-(2rabcos(phi))+(y+pphi)^2)
End Function
[/code]