' 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