' Positions of 1st 10 zeros in bj1(x)
data 3.83170597
data 7.01558668
data 10.17346814
data 13.32369194
data 16.47063005
data 19.61585851
data 22.76008438
data 25.90367209
data 29.04682854
data 32.18967991
'
dim zero(5)
for j=1 to 5
read zero(j)
next j
'
open "radial_grad.txt" for output as #1
max_x = zero(5)
np = 10
for j=0 to 4
for p=0 to np-1
x = zero(j)*(1-p/np) + zero(j+1)*p/np
if (x=0) then airy=1 : else airy = (2 * bj1(x) / x)^2
airy = airy^0.4 : ' Gamma function (increases brightness of low values)
x0 = x/max_x
?#1,"<stop offset=";chr$(34);str$(x0);chr$(34);
?#1," style=";chr$(34);"stop-color:#";
h$ = hex$(floor(255*abs(airy)+0.5),2)
?#1,h$;h$;h$;chr$(34);"/>"
next p
next j
?#1,"<stop offset=";chr$(34);"1";chr$(34);" style=";
?#1,chr$(34);"stop-color:#000000";chr$(34);"/>"
close#1
end
'
' BesselJ1 function based on Javascript source code at
' http://www.mhtl.uwaterloo.ca/old/courses/me3532/js/bessel.html
sub bj1(x)
ax = abs(x)
if (ax < 8.0) then
y=x*x
tmp = (-2972611.439+y*(15704.48260+y*(-30.16036606)))
ans1 = x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*tmp)))
ans2 = 144725228442.0+y*(2300535178.0+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0))))
ans = ans1/ans2
else
z=8.0/ax
y=z*z
xx=ax-2.356194491
ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4+y*(0.2457520174e-5+y*(-0.240337019e-6))))
ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+y*(-0.88228987e-6+y*0.105787412e-6)))
ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2)
if (x < 0.0) then ans = -ans
endif
bj1 = ans
end sub

open "linear_grad.txt" for output as #1
for x = 0 to 20
xx = ((20-x)/20)^2 : ' Create more stop colors in low intensity region
y = (xx)^0.4 : ' Gamma function
c$ = hex$(floor(255*y+0.5),2)
c$ = c$+c$+c$
print #1,"<stop offset=";chr$(34);str$(1-xx);chr$(34);
print #1," style=";chr$(34);" stop-color:#";c$;chr$(34);"/>"
next x
close #1
end