« previous next »
Pages: [1]   Go Down

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

 
thing19081008.png

 
2222.png

 
fdsafdsfdsfdsf.png

Logged
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!!
Logged
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!
Logged

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
Logged

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.
Logged
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

 
Logged
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!
Logged
Pages: [1]   Go Up
« previous next »