Author Topic: Domain Coloring  (Read 4953 times)

0 Members and 1 Guest are viewing this topic.

Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 1091
  • he lives
    • View Profile
Domain Coloring
« on: February 14, 2022, 01:54:11 pm »
Code: QB64: [Select]
  1. 'A = render anti-aliased image
  2. 'B = render standard image
  3.  
  4. 'MouseButton1 = Re-center
  5. 'MouseButton3 = Zoom x2 (in)
  6. 'MouseButton2 = Zoom /2 (out)
  7.  
  8. 'Leftarrow  = Next example plot
  9. 'Rightarrow = Previous example plot
  10. 'Uparrow    = Increase iterations (where applicable)
  11. 'DownArrow  = Decrease iterations (where applicable)
  12.  
  13. 'Num 1 = Shading option (HSV vs RGB vs Average)
  14. 'Num 2 = Shading scheme (1, 2, 3)
  15. 'Num 3 = Shade near origin
  16. 'Num 4 = Shade axes
  17. 'Num 5 = Shade integers
  18. 'Num 6 = Shade contours
  19. 'Num 7 = Greyscale
  20. 'Num 8 = Stencil
  21. 'Num 9 = Invert
  22. 'Num 0 = Bolden
  23.  
  24. _TITLE "Domain Coloring"
  25.  
  26. 'SCREEN _NEWIMAGE(250, 250, 32)
  27. SCREEN _NEWIMAGE(800, 800, 32)
  28. 'SCREEN _NEWIMAGE(1024, 768, 32)
  29. 'SCREEN _NEWIMAGE(1920, 1080, 32)
  30.  
  31. pi = 4 * ATN(1)
  32.  
  33. DIM AS DOUBLE zoom, xshift, yshift, re, im
  34. DIM AS STRING ExhibitName
  35. DIM SHARED PlotOption(0 TO 10) AS DOUBLE
  36.  
  37. PlotOption(0) = 10
  38. PlotOption(1) = .5
  39. PlotOption(2) = 1
  40. PlotOption(3) = 1
  41. PlotOption(4) = 1
  42. PlotOption(5) = 1
  43. PlotOption(6) = 1
  44. PlotOption(7) = -1
  45. PlotOption(8) = -1
  46. PlotOption(9) = -1
  47. PlotOption(10) = -1
  48.  
  49. ' Selector
  50. q = 1
  51.  
  52. Initialize q, ExhibitName, zoom, xshift, yshift
  53. CALL DrawPlot(ExhibitName, zoom, xshift, yshift, -1)
  54.  
  55. ' User interaction
  56. DIM ReDraw AS INTEGER
  57.         re = xshift + (_MOUSEX - _WIDTH / 2) / zoom
  58.         im = yshift - (_MOUSEY - _HEIGHT / 2) / zoom
  59.         Calculate re, im, ExhibitName
  60.         'LOCATE 5, 1: PRINT "re: "; INT(re + .5); "      "
  61.         'LOCATE 6, 1: PRINT "im: "; INT(im + .5); "      "
  62.         IF (_MOUSEBUTTON(1)) THEN
  63.             xshift = xshift + (_MOUSEX - _WIDTH / 2) / zoom
  64.             yshift = yshift - (_MOUSEY - _HEIGHT / 2) / zoom
  65.             ReDraw = -1
  66.             DO WHILE _MOUSEINPUT: LOOP
  67.         END IF
  68.         IF (_MOUSEBUTTON(2)) THEN
  69.             zoom = zoom * 1 / 2
  70.             ReDraw = -1
  71.             DO WHILE _MOUSEINPUT: LOOP
  72.         END IF
  73.         IF (_MOUSEBUTTON(3)) THEN
  74.             zoom = zoom * 2
  75.             ReDraw = -1
  76.             DO WHILE _MOUSEINPUT: LOOP
  77.         END IF
  78.     LOOP
  79.     kh = _KEYHIT
  80.     _KEYCLEAR
  81.     IF ((kh = ASC("a")) OR (kh = ASC("A"))) THEN
  82.         ReDraw = 1
  83.     END IF
  84.     IF ((kh = ASC("b")) OR (kh = ASC("B"))) THEN
  85.         ReDraw = -1
  86.     END IF
  87.     IF (kh = ASC("1")) THEN
  88.         PlotOption(1) = PlotOption(1) + .5
  89.         IF ((PlotOption(1) > 1)) THEN PlotOption(1) = 0
  90.         ReDraw = -1
  91.     END IF
  92.     IF (kh = ASC("2")) THEN
  93.         PlotOption(2) = PlotOption(2) + 1
  94.         IF (PlotOption(2) > 4) THEN PlotOption(2) = 1
  95.         ReDraw = -1
  96.     END IF
  97.     IF (kh = ASC("3")) THEN PlotOption(3) = -PlotOption(3): ReDraw = -1
  98.     IF (kh = ASC("4")) THEN PlotOption(4) = -PlotOption(4): ReDraw = -1
  99.     IF (kh = ASC("5")) THEN PlotOption(5) = -PlotOption(5): ReDraw = -1
  100.     IF (kh = ASC("6")) THEN PlotOption(6) = -PlotOption(6): ReDraw = -1
  101.     IF (kh = ASC("7")) THEN PlotOption(7) = -PlotOption(7): ReDraw = -1
  102.     IF (kh = ASC("8")) THEN PlotOption(8) = -PlotOption(8): ReDraw = -1
  103.     IF (kh = ASC("9")) THEN PlotOption(9) = -PlotOption(9): ReDraw = -1
  104.     IF (kh = ASC("0")) THEN PlotOption(10) = -PlotOption(10): ReDraw = -1
  105.     IF ((kh = ASC(" ")) OR (kh = 19712)) THEN
  106.         q = q + 1
  107.         Initialize q, ExhibitName, zoom, xshift, yshift
  108.         ReDraw = -1
  109.     END IF
  110.     IF (kh = 19200) THEN
  111.         q = q - 1
  112.         IF (q = 0) THEN q = 1
  113.         Initialize q, ExhibitName, zoom, xshift, yshift
  114.         ReDraw = -1
  115.     END IF
  116.     IF (kh = 18432) THEN
  117.         PlotOption(0) = PlotOption(0) + 1
  118.         ReDraw = -1
  119.     END IF
  120.     IF (kh = 20480) THEN
  121.         PlotOption(0) = PlotOption(0) - 1
  122.         IF (PlotOption(0) < 0) THEN PlotOption(0) = 0
  123.         ReDraw = -1
  124.     END IF
  125.     IF (ReDraw <> 0) THEN
  126.         CALL DrawPlot(ExhibitName, zoom, xshift, yshift, ReDraw)
  127.         ReDraw = 0
  128.     END IF
  129.  
  130.  
  131.  
  132. SUB Initialize (x AS LONG, ExhibitName AS STRING, zoom AS DOUBLE, xshift AS DOUBLE, yshift AS DOUBLE)
  133.     xshift = 0
  134.     yshift = 0
  135.     zoom = 50 * 2
  136.     SELECT CASE x
  137.         CASE 1
  138.             ExhibitName = "Vanilla"
  139.         CASE 2
  140.             ExhibitName = "Monomial"
  141.         CASE 3
  142.             ExhibitName = "Pole"
  143.         CASE 4
  144.             ExhibitName = "Shifts"
  145.         CASE 5
  146.             ExhibitName = "Geometric Series"
  147.             PlotOption(0) = 5
  148.         CASE 6
  149.             ExhibitName = "Taylor"
  150.             PlotOption(0) = 10
  151.         CASE 7
  152.             ExhibitName = "Exponential"
  153.         CASE 8
  154.             ExhibitName = "Logarithm"
  155.         CASE 9
  156.             ExhibitName = "Sqrt"
  157.         CASE 10
  158.             ExhibitName = "Branch"
  159.         CASE 11
  160.             ExhibitName = "Cosine"
  161.         CASE 12
  162.             ExhibitName = "Sine"
  163.         CASE 13
  164.             ExhibitName = "Gamma"
  165.         CASE 14
  166.             ExhibitName = "Condenser"
  167.         CASE 15
  168.             ExhibitName = "Inductor"
  169.         CASE 16
  170.             ExhibitName = "Mandelbrot"
  171.             xshift = -.5 - .25
  172.             PlotOption(0) = 1
  173.         CASE 17
  174.             ExhibitName = "Julia"
  175.             PlotOption(0) = 8
  176.         CASE 18
  177.             ExhibitName = "CIF1"
  178.             PlotOption(0) = 10
  179.         CASE 19
  180.             ExhibitName = "CIF2"
  181.             PlotOption(0) = 50
  182.         CASE 20
  183.             ExhibitName = "CIF3"
  184.         CASE 21
  185.             ExhibitName = "Canonical Logarithm"
  186.         CASE ELSE
  187.     END SELECT
  188.  
  189. SUB Calculate (x AS DOUBLE, y AS DOUBLE, a AS STRING)
  190.     DIM AS INTEGER j, m, n
  191.     DIM AS DOUBLE re, im, u, v, u0, v0, uu, vv, fu, fv, xx, yy, p, q, t
  192.  
  193.     SELECT CASE a
  194.         CASE "Vanilla"
  195.             re = x
  196.             im = y
  197.  
  198.         CASE "Monomial"
  199.             're = x ^ 2 - y ^ 2
  200.             'im = 2 * x * y
  201.             'cexp re, im, x, y, 2, 0
  202.             cexp re, im, x, y, 3, 0
  203.  
  204.         CASE "Pole"
  205.             'cexp re, im, x, y, -1, 0
  206.             cexp re, im, x, y, -2, 0
  207.  
  208.         CASE "Shifts"
  209.             'cdiv u, v, x + 1, y, x - 1, y
  210.             cdiv u, v, x, y, x - 1, y
  211.             re = u
  212.             im = v
  213.  
  214.         CASE "Geometric Series"
  215.             p = 0
  216.             q = 0
  217.             u = 1
  218.             v = 0
  219.             cadd u0, v0, p, q, u, v
  220.             FOR m = 1 TO PlotOption(0)
  221.                 cexp uu, vv, x, y, m, 0
  222.                 cadd p, q, u0, v0, uu, vv
  223.                 u0 = p
  224.                 v0 = q
  225.             NEXT
  226.             re = u0
  227.             im = v0
  228.  
  229.         CASE "Taylor"
  230.             p = 0
  231.             q = 0
  232.             u = 1
  233.             v = 0
  234.             cadd u0, v0, p, q, u, v
  235.             FOR m = 1 TO PlotOption(0)
  236.                 cexp uu, vv, x, y, m, 0
  237.                 cdiv u, v, uu, vv, facto&(m), 0
  238.                 cadd p, q, u0, v0, u, v
  239.                 u0 = p
  240.                 v0 = q
  241.             NEXT
  242.             re = u0
  243.             im = v0
  244.  
  245.         CASE "Exponential"
  246.             cexp re, im, EXP(1), 0, x, y
  247.  
  248.         CASE "Logarithm"
  249.             clog re, im, x, y
  250.  
  251.         CASE "Sqrt"
  252.             cexp re, im, x, y, 1 / 2, 0
  253.  
  254.         CASE "Branch"
  255.             cdiv p, q, x + 2, y + 1, x - 2, y - 1
  256.             clog re, im, p, q
  257.  
  258.         CASE "Cosine"
  259.             cosz re, im, x, y
  260.  
  261.         CASE "Sine"
  262.             sinz re, im, x, y
  263.  
  264.         CASE "Gamma"
  265.             ' Is this right?
  266.             uu = x - y
  267.             vv = x + y
  268.             cgamma uu, vv, x, y
  269.             re = uu
  270.             im = vv
  271.  
  272.         CASE "Condenser"
  273.             re = 1
  274.             im = 0
  275.             FOR j = 0 TO 7
  276.                 cmul u, v, 1, 0, -0.1 * (j - 3.5), 0
  277.                 cmul re, im, re, im, x - u, y - v - 0.1
  278.                 cdiv re, im, re, im, x - u, y - v + 0.1
  279.             NEXT
  280.  
  281.         CASE "Inductor"
  282.             re = 1
  283.             im = 0
  284.             FOR j = 0 TO 5
  285.                 cmul u, v, 1, 0, 0, -0.18 * (j - 2.5)
  286.                 cmul re, im, re, im, x - u - 0.2, y - v
  287.                 cdiv re, im, re, im, x - u + 0.2, y - v + 0.1
  288.             NEXT
  289.  
  290.         CASE "Mandelbrot"
  291.             u = x
  292.             v = y
  293.             FOR m = 0 TO PlotOption(0) '2 '5 '00
  294.                 u0 = u
  295.                 v0 = v
  296.                 u = u0 ^ 2 - v0 ^ 2 + x
  297.                 v = 2 * u0 * v0 + y
  298.             NEXT
  299.             re = u
  300.             im = v
  301.  
  302.         CASE "Julia"
  303.             u = x
  304.             v = y
  305.             FOR m = 0 TO PlotOption(0)
  306.                 u0 = u
  307.                 v0 = v
  308.                 u = u0 ^ 2 - v0 ^ 2 + 0.35
  309.                 v = 2 * u0 * v0 + 0 '0.5
  310.             NEXT
  311.             re = u
  312.             im = v
  313.  
  314.         CASE "CIF1"
  315.             ' Riemann sum resolution
  316.             n = PlotOption(0) '25
  317.             xx = 0
  318.             yy = 0
  319.  
  320.             FOR j = 0 TO n - 1
  321.  
  322.                 ' Integration contour
  323.                 ' 1:
  324.                 'u = 3 * (2 * COS(j * 2 * pi / n))
  325.                 'v = 3 * (2 * SIN(j * 2 * pi / n))
  326.                 ' 2:
  327.                 u = 4 * (2 * COS(j * 2 * pi / n))
  328.                 v = 4 * (2 * SIN(j * 2 * pi / n) * COS(j * 2 * pi / n))
  329.  
  330.                 ' f(z)
  331.                 fu = u - v 'u ^ 2 - v ^ 2
  332.                 fv = u + v '2 * u * v
  333.  
  334.                 ' f(z) / (z - z_0)
  335.                 cdiv uu, vv, fu, fv, u - x, v - y
  336.  
  337.                 ' z'(t) = derivative of integration contour
  338.                 ' 1:
  339.                 'cmul re, im, uu, vv, 3 * (-2 * SIN(j * 2 * pi / n)), 3 * (2 * COS(j * 2 * pi / n))
  340.                 ' 2:
  341.                 cmul re, im, uu, vv, 4 * (-2 * SIN(j * 2 * pi / n)), 4 * (2 * COS(j * 4 * pi / n))
  342.  
  343.                 ' Integral height calculation
  344.                 IF ((j = 0) OR (j = n - 1)) THEN
  345.                     xx = xx + 0.5 * re
  346.                     yy = yy + 0.5 * im
  347.                 ELSE
  348.                     xx = xx + re
  349.                     yy = yy + im
  350.                 END IF
  351.             NEXT
  352.  
  353.             ' Integral base calculation
  354.             xx = xx * 2 * pi / n
  355.             yy = yy * 2 * pi / n
  356.  
  357.             ' Rescale by 1/(2*pi*i)
  358.             cmul re, im, xx, yy, 0, -1 / (2 * pi)
  359.  
  360.         CASE "CIF2"
  361.             ' Riemann sum resolution
  362.             n = PlotOption(0) '10
  363.             xx = 0
  364.             yy = 0
  365.  
  366.             DIM xr AS DOUBLE
  367.             DIM xi AS DOUBLE
  368.             xr = 2.0
  369.             xi = 0.0
  370.  
  371.             FOR j = 0 TO n - 1
  372.  
  373.                 ' Integration contour
  374.                 u = 3 * 3 * COS(j * 2 * pi / n)
  375.                 v = 3 * 2 * SIN(j * 2 * pi / n)
  376.  
  377.                 ' f(z)
  378.                 sinz fu, fv, u, v
  379.  
  380.                 ' f(z) / (z - z_0)^(xr + i xi)
  381.                 p = u - x
  382.                 q = v - y
  383.                 cexp p, q, p, q, xr, xi
  384.                 cdiv uu, vv, fu, fv, p, q
  385.  
  386.                 ' z'(t) = derivative of integration contour
  387.                 cmul re, im, uu, vv, 3 * -3 * SIN(j * 2 * pi / n), 2 * COS(j * 2 * pi / n)
  388.  
  389.                 ' Integral height calculation
  390.                 IF ((j = 0) OR (j = n - 1)) THEN
  391.                     xx = xx + 0.5 * re
  392.                     yy = yy + 0.5 * im
  393.                 ELSE
  394.                     xx = xx + re
  395.                     yy = yy + im
  396.                 END IF
  397.             NEXT
  398.  
  399.             ' Integral base calculation
  400.             xx = xx * 2 * pi / n
  401.             yy = yy * 2 * pi / n
  402.  
  403.             ' rescale by Gamma(xr + i ri)/(2*pi*i)
  404.             cgamma u, v, xr + 1, xi
  405.             cmul re, im, xx, yy, u, v
  406.             cmul re, im, re, im, 0, -1 / (2 * pi)
  407.  
  408.         CASE "CIF3"
  409.             xx = 0
  410.             yy = 0
  411.             FOR t = 0 TO 3 STEP 0.1
  412.                 uu = t * COS(5 * t)
  413.                 vv = t * SIN(5 * t)
  414.                 cexp p, q, EXP(1), 0, uu, vv
  415.                 cdiv u0, v0, p, q, uu - x, vv - y
  416.                 cmul p, q, u0, v0, COS(5 * t) - 5 * t * SIN(5 * t), SIN(5 * t) + 5 * t * COS(5 * t)
  417.                 IF ((t = 0) OR (t = 3)) THEN
  418.                     xx = xx + 0.5 * p
  419.                     yy = yy + 0.5 * q
  420.                 ELSE
  421.                     xx = xx + 1 * p
  422.                     yy = yy + 1 * q
  423.                 END IF
  424.             NEXT
  425.             xx = xx * 0.01
  426.             yy = yy * 0.01
  427.             cmul re, im, xx, yy, 0, -1 / (2 * pi)
  428.  
  429.         CASE "Canonical Logarithm"
  430.             DIM AS DOUBLE a0, k0, sx0, sy0, tx0, ty0
  431.             a0 = .5
  432.             k0 = 0.09
  433.             sx0 = 2
  434.             sy0 = 1
  435.             tx0 = a0 * EXP(k0 * 16) * COS(16) + sx0
  436.             ty0 = a0 * EXP(k0 * 16) * SIN(16) + sy0
  437.             re = 0
  438.             im = 0
  439.             FOR t = 0 TO 16 STEP .1
  440.                 p = a0 * EXP(k0 * t) * COS(t) + sx0 - tx0
  441.                 q = a0 * EXP(k0 * t) * SIN(t) + sy0 - ty0
  442.                 cdiv uu, vv, 1, 0, p - x, q - y
  443.                 cmul xx, yy, uu, vv, a0 * EXP(k0 * t) * (k0 * COS(t) - SIN(t)), a0 * EXP(k0 * t) * (k0 * SIN(t) + COS(t))
  444.                 re = re + xx
  445.                 im = im + yy
  446.             NEXT
  447.             FOR t = 16 - .1 TO 0 STEP -.1
  448.                 p = -a0 * EXP(k0 * t) * COS(t) - sx0 + tx0
  449.                 q = -a0 * EXP(k0 * t) * SIN(t) - sy0 + ty0
  450.                 cdiv uu, vv, 1, 0, p - x, q - y
  451.                 cmul xx, yy, uu, vv, a0 * EXP(k0 * t) * (k0 * COS(t) - SIN(t)), a0 * EXP(k0 * t) * (k0 * SIN(t) + COS(t))
  452.                 re = re + xx
  453.                 im = im + yy
  454.             NEXT
  455.             re = re * 0.1
  456.             im = im * 0.1
  457.             cmul uu, vv, re, im, 0, -1 / (2 * pi)
  458.             re = uu
  459.             im = vv
  460.  
  461.         CASE ELSE
  462.     END SELECT
  463.  
  464.     x = re
  465.     y = im
  466.  
  467. SUB ShadePixel (red AS DOUBLE, grn AS DOUBLE, blu AS DOUBLE, alf AS DOUBLE, re AS DOUBLE, im AS DOUBLE, zoom AS DOUBLE)
  468.     DIM AS _UNSIGNED LONG c0, c1, c2
  469.     DIM AS DOUBLE r, a, h, s, k
  470.  
  471.     r = SQR(re * re + im * im)
  472.  
  473.     ' Color scheme 0
  474.     a = (pi + _ATAN2(im, -re)) / (2 * pi)
  475.     c1~& = hrgb~&(a, r)
  476.  
  477.     ' Color scheme 1
  478.     a = _ATAN2(im, -re) * .99999
  479.     h = 180 + a * 180 / pi
  480.     SELECT CASE PlotOption(2)
  481.         CASE 1
  482.             k = 50 * r
  483.         CASE 2
  484.             k = 50 * LOG(1 + r * zoom)
  485.         CASE 3
  486.             k = 50 * LOG(1 + r) * (1 + r)
  487.         CASE 4
  488.             k = 50 * r / zoom
  489.         CASE ELSE
  490.             'IF (r < 1) THEN
  491.             'ELSE
  492.             'END IF
  493.     END SELECT
  494.     s = 50 + ((k * 1) MOD 50)
  495.     c2~& = HSVtoRGB~&(h, s, 100)
  496.  
  497.     ' Weighted average of color schemes
  498.     c0~& = ShadeBlend(PlotOption(1), c1~&, c2~&)
  499.     c0~& = SetAlpha(c0~&, 255)
  500.  
  501.     'IF ((re > 0) AND (im > 0)) THEN
  502.  
  503.     ' Origin
  504.     IF (PlotOption(3) = 1) THEN c0~& = ShadeOrigin(c0~&, re, im, 0.075, 100)
  505.  
  506.     ' Axes
  507.     IF (PlotOption(4) = 1) THEN c0~& = ShadeAxes(c0~&, re, im, 0.075, 100)
  508.  
  509.     ' Integers
  510.     IF (PlotOption(5) = 1) THEN c0~& = ShadeIntegers(c0~&, re, im, 0.075, 100)
  511.  
  512.     'END IF
  513.  
  514.     ' Contours
  515.     IF (PlotOption(6) = 1) THEN c0~& = ShadeContours(c0~&, k, 50, 2, 100)
  516.  
  517.     ' Greyscale
  518.     IF (PlotOption(7) = 1) THEN c0~& = ShadeGreyscale(c0~&)
  519.  
  520.     ' Stencil
  521.     IF (PlotOption(8) = 1) THEN c0~& = ShadeStencil(c0~&, 100)
  522.  
  523.     ' Invert
  524.     IF (PlotOption(9) = 1) THEN c0~& = ShadeInvert(c0~&)
  525.  
  526.     ' Bolden
  527.     IF (PlotOption(0) = 1) THEN c0~& = ShadeBolden(c0~&, .8)
  528.  
  529.     red = _RED32(c0~&)
  530.     grn = _GREEN32(c0~&)
  531.     blu = _BLUE32(c0~&)
  532.     alf = _ALPHA32(c0~&)
  533.  
  534. SUB DrawPlot (TheExhibit AS STRING, zoom AS DOUBLE, xshift AS DOUBLE, yshift AS DOUBLE, AAtoggle AS INTEGER)
  535.     DIM AS DOUBLE j, k, x0, y0, re, im, red, grn, blu, alf, nr, ng, nb, na, d, f, jj, kk
  536.     DIM AS DOUBLE r(99), g(99), b(99), a(99), w(99)
  537.     DIM AS INTEGER ii
  538.  
  539.     CLS
  540.     _TITLE TheExhibit
  541.     FOR j = 0 TO _WIDTH
  542.         FOR k = 0 TO _HEIGHT
  543.  
  544.             IF (AAtoggle = -1) THEN
  545.                 x0 = (j) - _WIDTH / 2
  546.                 y0 = -(k) + _HEIGHT / 2
  547.                 re = x0 / zoom + xshift
  548.                 im = y0 / zoom + yshift
  549.                 Calculate re, im, TheExhibit
  550.                 ShadePixel red, grn, blu, alf, re, im, zoom
  551.  
  552.             ELSE
  553.                 f = 2
  554.                 d = .25
  555.                 ii = 0
  556.                 FOR jj = j - f * d TO j + f * d STEP d
  557.                     FOR kk = k - f * d TO k + f * d STEP d
  558.                         ii = ii + 1
  559.                         x0 = (jj - d) - _WIDTH / 2
  560.                         y0 = -(kk - d) + _HEIGHT / 2
  561.                         re = x0 / zoom + xshift
  562.                         im = y0 / zoom + yshift
  563.                         Calculate re, im, TheExhibit
  564.                         ShadePixel red, grn, blu, alf, re, im, zoom
  565.                         r(ii) = red
  566.                         g(ii) = grn
  567.                         b(ii) = blu
  568.                         a(ii) = alf
  569.                         w(ii) = EXP(-1 * ((j - jj) ^ 2 + (k - kk) ^ 2))
  570.                     NEXT
  571.                 NEXT
  572.                 red = 0
  573.                 grn = 0
  574.                 blu = 0
  575.                 alf = 0
  576.                 nr = 0
  577.                 ng = 0
  578.                 nb = 0
  579.                 na = 0
  580.                 FOR jj = 1 TO ii
  581.                     red = red + r(jj) * w(ii)
  582.                     grn = grn + g(jj) * w(ii)
  583.                     blu = blu + b(jj) * w(ii)
  584.                     alf = alf + a(jj) * w(ii)
  585.                     nr = nr + w(ii)
  586.                     ng = ng + w(ii)
  587.                     nb = nb + w(ii)
  588.                     na = na + w(ii)
  589.                 NEXT
  590.                 red = red / nr
  591.                 grn = grn / ng
  592.                 blu = blu / nb
  593.                 alf = alf / na
  594.             END IF
  595.  
  596.             CALL CPset(x0, y0, _RGBA(red, grn, blu, alf))
  597.  
  598.         NEXT
  599.     NEXT
  600.  
  601.     'COLOR _RGB32(255, 255, 255, 255)
  602.     'LOCATE 1, 1: PRINT TheExhibit
  603.     'LOCATE 2, 1: PRINT "zoom: "; zoom
  604.     'LOCATE 3, 1: PRINT "xshift: "; xshift
  605.     'LOCATE 4, 1: PRINT "yshift: "; yshift
  606.  
  607. FUNCTION ShadeBlend~& (f AS DOUBLE, x AS _UNSIGNED LONG, y AS _UNSIGNED LONG)
  608.     DIM AS DOUBLE red, grn, blu, alf
  609.     red = (1 - f) * (_RED32(x)) + f * (_RED32(y))
  610.     grn = (1 - f) * (_GREEN32(x)) + f * (_GREEN32(y))
  611.     blu = (1 - f) * (_BLUE32(x)) + f * (_BLUE32(y))
  612.     alf = (1 - f) * (_ALPHA32(x)) + f * (_ALPHA32(y))
  613.     ShadeBlend~& = _RGB32(red, grn, blu, alf)
  614.  
  615. FUNCTION SetAlpha~& (x AS _UNSIGNED LONG, y AS DOUBLE)
  616.     SetAlpha~& = _RGB32(_RED32(x), _GREEN32(x), _BLUE32(x), y)
  617.  
  618. FUNCTION ShadeContours~& (c AS _UNSIGNED LONG, k AS DOUBLE, x AS DOUBLE, y AS DOUBLE, a AS DOUBLE)
  619.     IF (((k MOD x) < y) OR ((k MOD x) > (x - y))) THEN c = SetAlpha(c, a)
  620.     ShadeContours~& = c
  621.  
  622. FUNCTION ShadeIntegers~& (c AS _UNSIGNED LONG, x AS DOUBLE, y AS DOUBLE, z AS DOUBLE, a AS DOUBLE)
  623.     IF (((x - INT(x)) < z) AND (INT(x) <> 0)) THEN c = SetAlpha(c, a)
  624.     IF (((y - INT(y)) < z) AND (INT(y) <> 0)) THEN c = SetAlpha(c, a)
  625.     IF (((INT(x) + 1 - x) < z) AND ((INT(x) + 1) <> 0)) THEN c = SetAlpha(c, a)
  626.     IF (((INT(y) + 1 - y) < z) AND ((INT(y) + 1) <> 0)) THEN c = SetAlpha(c, a)
  627.     ShadeIntegers~& = c
  628.  
  629. FUNCTION ShadeAxes~& (c AS _UNSIGNED LONG, x AS DOUBLE, y AS DOUBLE, z AS DOUBLE, a AS DOUBLE)
  630.     IF SQR(x * x + y * y) > 1 THEN
  631.         IF (((x - INT(x)) < z) AND (INT(x) = 0)) THEN c = SetAlpha(c, a)
  632.         IF (((y - INT(y)) < z) AND (INT(y) = 0)) THEN c = SetAlpha(c, a)
  633.         IF (((INT(x) + 1 - x) < z) AND ((INT(x) + 1) = 0)) THEN c = SetAlpha(c, a)
  634.         IF (((INT(y) + 1 - y) < z) AND ((INT(y) + 1) = 0)) THEN c = SetAlpha(c, a)
  635.     END IF
  636.     ShadeAxes~& = c
  637.  
  638. FUNCTION ShadeOrigin~& (c AS _UNSIGNED LONG, x AS DOUBLE, y AS DOUBLE, z AS DOUBLE, a AS DOUBLE)
  639.     IF SQR(x * x + y * y) <= 1 THEN
  640.         IF (((x - INT(x)) < z) AND (INT(x) = 0)) THEN c = SetAlpha(c, a)
  641.         IF (((y - INT(y)) < z) AND (INT(y) = 0)) THEN c = SetAlpha(c, a)
  642.         IF (((INT(x) + 1 - x) < z) AND ((INT(x) + 1) = 0)) THEN c = SetAlpha(c, a)
  643.         IF (((INT(y) + 1 - y) < z) AND ((INT(y) + 1) = 0)) THEN c = SetAlpha(c, a)
  644.     END IF
  645.     ShadeOrigin~& = c
  646.  
  647. FUNCTION ShadeGreyscale~& (x AS _UNSIGNED LONG)
  648.     DIM AS DOUBLE a
  649.     a = 255 - (1 / 3) * (_RED32(x) + _GREEN32(x) + _BLUE32(x))
  650.     ShadeGreyscale~& = _RGB32(a, a, a, _ALPHA32(x))
  651.  
  652. FUNCTION ShadeStencil~& (x AS _UNSIGNED LONG, t AS DOUBLE)
  653.     DIM AS DOUBLE red, grn, blu
  654.     IF (_ALPHA32(x) <> t) THEN
  655.         red = 0
  656.         grn = 0
  657.         blu = 0
  658.     ELSE
  659.         red = _RED32(x)
  660.         grn = _GREEN32(x)
  661.         blu = _BLUE32(x)
  662.     END IF
  663.     ShadeStencil~& = _RGB32(red, grn, blu, 255)
  664.  
  665. FUNCTION ShadeInvert~& (x AS _UNSIGNED LONG)
  666.     ShadeInvert~& = _RGB32(255 - _RED32(x), 255 - _GREEN32(x), 255 - _BLUE32(x), _ALPHA32(x))
  667.  
  668. FUNCTION ShadeBolden~& (x AS _UNSIGNED LONG, t AS DOUBLE)
  669.     DIM AS DOUBLE red, grn, blu
  670.     red = _RED32(x)
  671.     grn = _GREEN32(x)
  672.     blu = _BLUE32(x)
  673.     IF ((red + grn + blu) < (3 * 255)) THEN
  674.         red = red * t
  675.         grn = grn * t
  676.         blu = blu * t
  677.     END IF
  678.     ShadeBolden~& = _RGB32(red, grn, blu, _ALPHA32(x))
  679.  
  680. FUNCTION HSVtoRGB~& (h AS DOUBLE, s AS DOUBLE, v AS DOUBLE)
  681.     DIM AS DOUBLE c, x, m, r, g, b
  682.     IF ((h > 360) OR (h < 0) OR (s > 100) OR (s < 0) OR (v > 100) OR (v < 0)) THEN
  683.         PRINT "Out of range:"; h; s; v
  684.     END IF
  685.     s = s / 100
  686.     v = v / 100
  687.     c = s * v
  688.     x = c * (1 - ABS(fmod(h / 60, 2) - 1))
  689.     m = v - c
  690.     IF ((h >= 0) AND (h < 60)) THEN
  691.         r = c
  692.         g = x
  693.         b = 0
  694.     END IF
  695.     IF ((h >= 60) AND (h < 120)) THEN
  696.         r = x
  697.         g = c
  698.         b = 0
  699.     END IF
  700.     IF ((h >= 120) AND (h < 180)) THEN
  701.         r = 0
  702.         g = c
  703.         b = x
  704.     END IF
  705.     IF ((h >= 180) AND (h < 240)) THEN
  706.         r = 0
  707.         g = x
  708.         b = c
  709.     END IF
  710.     IF ((h >= 240) AND (h < 300)) THEN
  711.         r = x
  712.         g = 0
  713.         b = c
  714.     END IF
  715.     IF ((h >= 300) AND (h < 360)) THEN
  716.         r = c
  717.         g = 0
  718.         b = x
  719.     END IF
  720.     r = (r + m) * 255
  721.     g = (g + m) * 255
  722.     b = (b + m) * 255
  723.     HSVtoRGB~& = _RGB32(r, g, b)
  724.  
  725. FUNCTION fmod## (numer AS DOUBLE, denom AS DOUBLE)
  726.     fmod## = numer - INT(numer / denom) * denom
  727.  
  728. ' a vince original, slightly modified
  729. FUNCTION hrgb~& (h AS DOUBLE, m AS DOUBLE)
  730.     DIM n AS LONG
  731.     DIM AS DOUBLE r, g, b, mm, p, rr, gg, bb
  732.     r = 0.5 - 0.5 * SIN(2 * pi * h - pi / 2)
  733.     g = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 - pi / 2)) * -(h < 0.66)
  734.     b = (0.5 + 0.5 * SIN(2 * pi * h * 1.5 + pi / 2)) * -(h > 0.33)
  735.     mm = 0 * (m * 500 MOD 500)
  736.     n = 16
  737.     p = ABS((h * n) - INT(h * n))
  738.     rr = 255 * r - 0.15 * mm - 35 * p
  739.     gg = 255 * g - 0.15 * mm - 35 * p
  740.     bb = 255 * b - 0.15 * mm - 35 * p
  741.     IF (rr < 0) THEN rr = 0
  742.     IF (gg < 0) THEN gg = 0
  743.     IF (bb < 0) THEN bb = 0
  744.     hrgb~& = _RGB(rr, gg, bb)
  745.  
  746. SUB CPset (x0 AS DOUBLE, y0 AS DOUBLE, shade AS _UNSIGNED LONG)
  747.     PSET (_WIDTH / 2 + x0, -y0 + _HEIGHT / 2), shade
  748.  
  749. 'SUB circlef (x AS LONG, y AS LONG, r AS LONG, c AS LONG)
  750. '    DIM AS LONG xx, yy, e
  751. '    xx = r
  752. '    yy = 0
  753. '    e = -r
  754. '    DO WHILE (yy < xx)
  755. '        IF (e <= 0) THEN
  756. '            yy = yy + 1
  757. '            LINE (x - xx, y + yy)-(x + xx, y + yy), c, BF
  758. '            LINE (x - xx, y - yy)-(x + xx, y - yy), c, BF
  759. '            e = e + 2 * yy
  760. '        ELSE
  761. '            LINE (x - yy, y - xx)-(x + yy, y - xx), c, BF
  762. '            LINE (x - yy, y + xx)-(x + yy, y + xx), c, BF
  763. '            xx = xx - 1
  764. '            e = e - 2 * xx
  765. '        END IF
  766. '    LOOP
  767. '    LINE (x - r, y)-(x + r, y), c, BF
  768. 'END SUB
  769.  
  770. SUB cadd (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE, aa AS DOUBLE, bb AS DOUBLE)
  771.     DIM AS DOUBLE x, y, a, b
  772.     x = xx
  773.     y = yy
  774.     a = aa
  775.     b = bb
  776.     u = x + a
  777.     v = y + b
  778.  
  779. SUB cmul (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE, aa AS DOUBLE, bb AS DOUBLE)
  780.     DIM AS DOUBLE x, y, a, b
  781.     x = xx
  782.     y = yy
  783.     a = aa
  784.     b = bb
  785.     u = x * a - y * b
  786.     v = x * b + y * a
  787.  
  788. SUB cdiv (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE, aa AS DOUBLE, bb AS DOUBLE)
  789.     DIM AS DOUBLE x, y, a, b, d
  790.     x = xx
  791.     y = yy
  792.     a = aa
  793.     b = bb
  794.     d = a * a + b * b
  795.     u = (x * a + y * b) / d
  796.     v = (y * a - x * b) / d
  797.  
  798. SUB cexp (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE, aa AS DOUBLE, bb AS DOUBLE)
  799.     DIM AS DOUBLE x, y, a, b, lnz, argz, mag, ang
  800.     x = xx
  801.     y = yy
  802.     a = aa
  803.     b = bb
  804.     lnz = x * x + y * y
  805.     IF (lnz = 0) THEN
  806.         u = 0
  807.         v = 0
  808.     ELSE
  809.         lnz = 0.5 * LOG(lnz)
  810.         argz = _ATAN2(y, x)
  811.         mag = EXP(a * lnz - b * argz)
  812.         ang = a * argz + b * lnz
  813.         u = mag * COS(ang)
  814.         v = mag * SIN(ang)
  815.     END IF
  816.  
  817. SUB clog (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE)
  818.     DIM AS DOUBLE x, y, lnz, argz
  819.     x = xx
  820.     y = yy
  821.     lnz = x * x + y * y
  822.     IF (lnz = 0) THEN
  823.         u = 0
  824.         v = 0
  825.     ELSE
  826.         lnz = 0.5 * LOG(lnz)
  827.         argz = _ATAN2(y, x)
  828.         u = lnz
  829.         v = argz
  830.     END IF
  831.  
  832. FUNCTION cosh## (x AS DOUBLE)
  833.     cosh## = 0.5## * (EXP(x) + EXP(-x))
  834.  
  835. FUNCTION sinh## (x AS DOUBLE)
  836.     sinh## = 0.5## * (EXP(x) - EXP(-x))
  837.  
  838. SUB sinz (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE)
  839.     DIM AS DOUBLE x, y
  840.     x = xx
  841.     y = yy
  842.     u = SIN(x) * cosh(y)
  843.     v = COS(x) * sinh(y)
  844.  
  845. SUB cosz (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE)
  846.     DIM AS DOUBLE x, y
  847.     x = xx
  848.     y = yy
  849.     u = COS(x) * cosh(y)
  850.     v = -SIN(x) * sinh(y)
  851.  
  852. FUNCTION rgamma## (x)
  853.     rgamma = SQR(2 * pi * x) * ((x / EXP(1)) ^ x)
  854.  
  855. SUB cgamma (u AS DOUBLE, v AS DOUBLE, xx AS DOUBLE, yy AS DOUBLE)
  856.     DIM AS DOUBLE x, y, uu, vv, p, q
  857.     x = xx
  858.     y = yy
  859.     IF x = 0 AND y = 0 THEN
  860.         u = 1
  861.         v = 0
  862.     ELSE
  863.         cexp uu, vv, x, y, x - 0.5, y
  864.         cexp p, q, EXP(1), 0, -x, -y
  865.         cmul u, v, uu, vv, p, q
  866.         u = u * SQR(2 * pi)
  867.         v = v * SQR(2 * pi)
  868.     END IF
  869.  
  870. FUNCTION facto& (x AS LONG)
  871.     IF (x = 1) THEN facto& = 1
  872.     IF (x > 1) THEN facto& = x * facto&(x - 1)
  873.  

 
thing19081008.png

 
2222.png

 
fdsafdsfdsfdsf.png

You're not done when it works, you're done when it's right.

Offline johnno56

  • Forum Resident
  • Posts: 1270
  • Live long and prosper.
    • View Profile
Re: Domain Coloring
« Reply #1 on: February 14, 2022, 02:45:20 pm »
Very cool, indeed!!
Logic is the beginning of wisdom.

Offline bplus

  • Global Moderator
  • Forum Resident
  • Posts: 8053
  • b = b + ...
    • View Profile
Re: Domain Coloring
« Reply #2 on: February 14, 2022, 03:00:58 pm »
+1 Nice!

Offline Dav

  • Forum Resident
  • Posts: 792
    • View Profile
Re: Domain Coloring
« Reply #3 on: February 14, 2022, 04:57:29 pm »
Really nice, @STxAxTIC!

- Dav

Offline STxAxTIC

  • Library Staff
  • Forum Resident
  • Posts: 1091
  • he lives
    • View Profile
Re: Domain Coloring
« Reply #4 on: February 14, 2022, 06:47:15 pm »
I should say, the whole reason this program exists is vince. This isn't the first domain coloring program out here, not even the first this week - and I even used his examples as a checkpoint to make sure this code was doing its job. He knows this though - we collaborated heavily. So definitely tilt the hats in vince's direction, all I did was overengineer his original idea.
You're not done when it works, you're done when it's right.

Offline Pete

  • Forum Resident
  • Posts: 2361
  • Cuz I sez so, varmint!
    • View Profile
Re: Domain Coloring
« Reply #5 on: February 25, 2022, 07:22:26 pm »
For the function...

FUNCTION facto& (x AS LONG)

if x = 1 then facto& = 1 else facto& = x * facto&(x - 1)

Anyway, don't you just hate it when we have to deal with an extra statement to handle the times zero conundrum?

Love the coding structure.

You might want to consider using it for a computer modern abstract art challenge. Who can create the most striking image in both shape and color. Just a thought. I have a few others, but they are not suitable for a family forum. 

Pete

 
Want to learn how to write code on cave walls? https://www.tapatalk.com/groups/qbasic/qbasic-f1/

Offline Phlashlite

  • Newbie
  • Posts: 50
    • View Profile
Re: Domain Coloring
« Reply #6 on: February 27, 2022, 04:20:52 am »
Very nice!