How to draw one sunset image using Mathematica?

(* ::Package:: *)
(* Procedural sunset simulation in Wolfram Language.

   Run from PowerShell:
   & "D:\Software\Wolfram Research\Mathematica\14.0\wolframscript.exe" -script .\sunset_simulation.wl
*)

SeedRandom[20260430];

width = 1280;
height = 720;
horizon = 0.43;              (* y coordinate, 0 = bottom, 1 = top *)
sunX = 0.66;
sunY = 0.49;
sunR = 0.088;

smoothstep[a_, b_, x_] := Module[{t = Clip[(x - a)/(b - a), {0, 1}]},
  t*t*(3 - 2*t)
];

mix[a_, b_, t_] := (1 - t) a + t b;
clampColor[c_] := Clip[c, {0, 1}];

skyColor[x_, y_] := Module[
  {
    t, base, d, glow, disk, cloudNoise1, cloudNoise2, clouds,
    horizonCol, midCol, topCol, cloudCol
  },
  t = Clip[(y - horizon)/(1 - horizon), {0, 1}];

  horizonCol = {1.00, 0.43, 0.16};
  midCol = {0.78, 0.18, 0.35};
  topCol = {0.07, 0.08, 0.22};
  cloudCol = {0.11, 0.07, 0.15};

  base = If[
    t < 0.45,
    mix[horizonCol, midCol, smoothstep[0, 0.45, t]],
    mix[midCol, topCol, smoothstep[0.35, 1, t]]
  ];

  d = Sqrt[((x - sunX)*width/height)^2 + (y - sunY)^2];
  glow = Exp[-(d/0.23)^2];
  disk = smoothstep[sunR + 0.006, sunR - 0.006, d];

  base = base + 0.36 glow {1.0, 0.42, 0.10};
  base = mix[base, {1.0, 0.77, 0.23}, disk];

  cloudNoise1 =
    0.52 + 0.22 Sin[28 x + 8 y] +
    0.18 Sin[55 x - 19 y + 1.2] +
    0.12 Sin[93 x + 7 y + 0.4];
  cloudNoise2 =
    0.49 + 0.25 Sin[35 x - 10 y + 2.1] +
    0.18 Sin[72 x + 27 y] +
    0.10 Sin[120 x - 13 y];

  clouds = Max[
    smoothstep[0.55, 0.78, cloudNoise1] Exp[-((y - 0.73)/0.060)^2],
    0.78 smoothstep[0.54, 0.76, cloudNoise2] Exp[-((y - 0.61)/0.044)^2]
  ];

  base = mix[base, cloudCol, 0.48 clouds];
  base = base + 0.13 clouds glow {1.0, 0.52, 0.18};

  clampColor[base]
];

waterColor[x_, y_] := Module[
  {
    depth, base, wave1, wave2, waveShade, spread, core, sparkle,
    reflectedSun, vignette
  },
  depth = Clip[(horizon - y)/horizon, {0, 1}];
  base = mix[{0.96, 0.34, 0.13}, {0.02, 0.07, 0.15}, smoothstep[0, 1, depth]];

  wave1 = 0.5 + 0.5 Sin[420 y + 16 Sin[14 x] + 35 depth x];
  wave2 = 0.5 + 0.5 Sin[860 y + 45 x + 8 Sin[31 x]];
  waveShade = 0.78 + 0.16 wave1 + 0.08 wave2;

  spread = 0.030 + 0.34 depth;
  core = Exp[-((x - sunX)/spread)^2] Exp[-2.65 depth];
  sparkle = Max[0, Sin[650 y + 42 Sin[23 x] + 11 Sin[70 x]]]^8;
  reflectedSun = core (0.26 + 1.65 sparkle);

  base = waveShade base + reflectedSun {1.0, 0.58, 0.13};
  vignette = 1 - 0.20 Abs[x - 0.5] - 0.10 depth;

  clampColor[vignette base]
];

pixelColor[col_, row_] := Module[{x, y, c},
  x = (col - 0.5)/width;
  y = 1 - (row - 0.5)/height;
  c = If[y >= horizon, skyColor[x, y], waterColor[x, y]];
  clampColor[c]
];

Print["Rendering procedural sky and water..."];
baseImage = Image[
  Table[pixelColor[col, row], {row, 1, height}, {col, 1, width}],
  "Real"
];

horizonPx = height horizon;

mountainTop[x_] :=
  horizonPx + 18 +
  28 Sin[2 Pi (1.15 x/width) + 0.6] +
  13 Sin[2 Pi (3.70 x/width) + 1.2] +
  8 Cos[2 Pi (7.10 x/width)];

landTop[x_] :=
  48 + 16 Sin[2 Pi (1.8 x/width) + 0.2] +
  9 Sin[2 Pi (5.3 x/width) + 1.7];

mountainPoints = Table[{x, mountainTop[x]}, {x, 0, width, 16}];
mountainPolygon = Join[
  {{0, horizonPx - 18}, {width, horizonPx - 18}},
  Reverse[mountainPoints]
];

landPoints = Table[{x, landTop[x]}, {x, 0, width, 16}];
landPolygon = Join[{{0, 0}, {width, 0}}, Reverse[landPoints]];

bird[x_, y_, s_] := {
  AbsoluteThickness[2.1],
  CapForm["Round"],
  BezierCurve[{{x - 2 s, y}, {x - s, y + 0.75 s}, {x, y}}],
  BezierCurve[{{x, y}, {x + s, y + 0.75 s}, {x + 2 s, y}}]
};

palmLeaves[x_, y_, s_] := Table[
  {
    AbsoluteThickness[9],
    CapForm["Round"],
    BezierCurve[{
      {x, y},
      {x + s 70 Cos[a], y + s 34 Sin[a] + 15},
      {x + s 130 Cos[a], y + s 52 Sin[a]}
    }]
  },
  {a, {0.25, 0.65, 1.05, 1.45, 2.05, 2.65, 3.15, 3.65}}
];

overlayGraphic = Graphics[
  {
    RGBColor[0.035, 0.026, 0.050],
    EdgeForm[None],
    Polygon[mountainPolygon],

    RGBColor[0.020, 0.018, 0.032],
    AbsoluteThickness[3],
    Line[{{0, horizonPx - 18}, {width, horizonPx - 18}}],

    RGBColor[0.025, 0.020, 0.035],
    Polygon[landPolygon],

    RGBColor[0.018, 0.015, 0.026],
    AbsoluteThickness[19],
    CapForm["Round"],
    BezierCurve[{{142, 46}, {154, 126}, {171, 246}}],
    AbsoluteThickness[13],
    BezierCurve[{{175, 54}, {187, 128}, {201, 226}}],
    palmLeaves[171, 246, 1.0],
    palmLeaves[201, 226, 0.72],

    RGBColor[0.018, 0.015, 0.026],
    AbsoluteThickness[10],
    Line[{{870, 56}, {1130, 92}}],
    AbsoluteThickness[7],
    Line[{{940, 49}, {952, 92}}],
    Line[{{1020, 48}, {1033, 82}}],
    Line[{{1110, 46}, {1119, 68}}],

    RGBColor[0.055, 0.040, 0.065],
    bird[815, 512, 12],
    bird[860, 535, 8],
    bird[902, 501, 10],
    bird[760, 548, 7]
  },
  PlotRange -> {{0, width}, {0, height}},
  ImageSize -> {width, height},
  AspectRatio -> height/width,
  Background -> None
];

Print["Compositing silhouettes..."];
baseGraphic = Graphics[
  {Inset[baseImage, {width/2, height/2}, Center, {width, height}]},
  PlotRange -> {{0, width}, {0, height}},
  ImageSize -> {width, height},
  AspectRatio -> height/width,
  Background -> Black
];

finalGraphic = Show[
  baseGraphic,
  overlayGraphic,
  PlotRange -> {{0, width}, {0, height}},
  ImageSize -> {width, height},
  AspectRatio -> height/width
];

finalImage = Rasterize[
  finalGraphic,
  "Image",
  ImageSize -> {width, height},
  RasterSize -> {width, height}
];

outputPath = FileNameJoin[{Directory[], "sunset_mathematica.png"}];
Export[outputPath, finalImage, "PNG"];

Print["Exported: " <> outputPath];

The sunset image is generated procedurally: no photo is used. The image plane is treated as normalized coordinates where $x \in [0,1]$ runs left to right and $y \in [0,1]$ runs bottom to top. A horizon value $h=0.43$ separates sky from water, so the pixel color is chosen by the piecewise rule

$$ C(x,y)= \begin{cases} C_{\text{sky}}(x,y), & y \ge h,\\ C_{\text{water}}(x,y), & y < h. \end{cases} $$

The sky is mainly a vertical color field. Near the horizon it is orange, then it blends through red and purple into dark blue near the top. The blend uses linear interpolation,

$$ \operatorname{mix}(a,b,t)=(1-t)a+tb, $$

where $a$ and $b$ are RGB colors and $t$ is a normalized height above the horizon. Soft transitions, such as the sun edge and cloud masks, use the common graphics function

$$ \operatorname{smoothstep}(t)=t^2(3-2t), $$

which avoids hard edges by easing smoothly from 0 to 1.

The sun is modeled as a disk plus a radial glow. If the sun center is $(x_s,y_s)$, the distance from a pixel to the sun is approximately

$$ d=\sqrt{\left((x-x_s)\frac{W}{H}\right)^2+(y-y_s)^2}, $$

where $W/H$ corrects for the wide image aspect ratio. The warm glow is then generated with an exponential falloff such as $g=e^{-(d/r)^2}$, so pixels close to the sun become bright yellow-orange while distant pixels receive only a weak tint.

The water uses the same idea, but with wave and reflection terms. A depth value $q=(h-y)/h$ darkens the water as it moves toward the viewer. Horizontal ripples are produced by sine waves such as $\sin(420y+16\sin(14x))$. The sun reflection is strongest below the sun and fades sideways and downward:

$$ R(x,y)=\exp\left(-\frac{(x-x_s)^2}{s(q)^2}\right)e^{-kq}. $$

Extra bright sparkles are made by taking only positive sine-wave peaks and raising them to a high power, which leaves thin highlights on the water surface.

After the pixel-based sky and water are rendered, vector shapes are composited on top. Mountains, land, palm trees, birds, and a pier are drawn with Polygon, Line, and BezierCurve. Because these shapes use very dark colors, they read as silhouettes against the bright sunset. The full pipeline is therefore: compute a color for every pixel, convert the color table into an image, draw vector silhouettes over it, rasterize the final composition, and export sunset_mathematica.png.