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/(pi*2*r))

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+sy*dw*

rb=r+tydw

L=Msnow(n,pitch,ra,rb,y,MaxErr)

'Lc=Msnow(n,pitch,ra,rb,yc,MaxErr)

rra=ra*ra*

rrb=rbrb

p=pitch/2/pi

LeftTermC=1e-7*sqrt(rra+p*p)*sqrt(rrb+p*p)/sqrt((1+p*p/rra)*(1+p*p/rrb))*

LeftTermU=1e-7sqrt(rra*rrb)/sqrt((1+p*p/rra)*(1+p*p/rrb))

cIterTerm=LeftTermC*L*

uIterTerm=LeftTermUL

montesumC=montesumC+cIterTerm

montesumU=montesumU+uIterTerm

cSqSum=cSqSum+cIterTerm*cIterTerm*

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=2*pi*N

nstr=""

rra=ra*ra*

rrb=rbrb

rab=ra*rb*

LeftTermC=1e-7sqrt(rra+p*p)**sqrt(rrb+p*p)/sqrt((1+pp/rra)*(1+p*p/rrb))

LeftTermU=1e-7*sqrt(rra*rrb)/sqrt((1+p*p/rra)*(1+p*p/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/2*dx*

While CurrentErr>MaxErr or m<512

m=2m

dx=dx/2

Sum=0

for i=1 to m step 2

Sum=Sum+HlxIntgrndRad(1,-i*dx-a,Ntp,p,rra,rrb,rab,y)+HlxIntgrndRad(-1,i*dx+a,Ntp,p,rra,rrb,rab,y)

next

Integral=(4*Sum+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+sn*phi)*(cos(phi)+p*p/rab)/sqrt(rra+rrb-(2*rab*cos(phi))+(y+p*phi)^2)

End Function

[/code]