2013/10/19

C言語とgnuplotでマンデルブロ集合を表示する

gnuplotを使うと結構簡単にマンデルブロ集合を表示することができます。OSはdebianを使用しています。

リスト1 mandelblot.c
/*
  マンデルブロ集合計算プログラム

 */
#include <stdio.h>
#include <math.h>

#define C0r     -0.743  // 計算する場所の中心座標(実数部)
#define C0i      0.1145 // 計算する場所の中心座標(虚数部)
#define VS       0.003  // 計算する場所の中心座標からの範囲(±VS)

#define NMAX     20000  // 計算の繰り返し上限
#define STEP     800.0  // 計算する刻み


double mandelbrot(double a, double b){
  double x = 0.0;
  double y = 0.0;
  double x1, y1;

  int n;

  for (n = 1; n <= NMAX; n++) {
    x1 = x * x - y * y + a;
    y1 = 2.0 * x * y + b;
    if ( x1 * x1 + y1 * y1 > 4.0) return log(n); // 発散
    x = x1;
    y = y1;
  }
  return 0; // 計算の繰り返し上限到達
}


int main() {
  double a, b;

  for (a = C0r-VS; a < C0r+VS; a += 2.0*VS/STEP) {
    for (b = C0i-VS; b < C0i+VS; b += 2.0*VS/STEP) {
      printf("%1.14e %1.14e %1.14e\n", a, b, mandelbrot(a, b));
    }
    printf("\n"); // これがないとgnuplotでエラーが出る
  }
  return 0;
}

リスト2 png.gnuplot (gnuplot実行スクリプト)
set term png size 900, 900
set output "mandelbrot.png"
set grid
set pm3d map
set size square
set palette defined (0 "#000000", 2 "#c00000", 7 "#ffff00", 9 "#ffffff")
splot "data.txt"

リスト3 コンパイル&実行
# cc -o mandelbrot -lm -O4 -Wall mandelbrot.c
# ./mandelbrot > data.txt
# gnuplot png.gnuplot

図1 実行結果(mandelbrot.png)
  • 計算する場所はリスト1のdefine文の数値を変えることにより可能です。
  • 色の変化を分かりやすくするため、Logスケールで表示しています。
C言語のほかにPythonとR言語で記述した記事もあります。詳しくは マンデルブロ集合 をご覧ください。

0 件のコメント:

コメントを投稿

【備忘録】時系列データの編集方法(R言語, tidyverse)

TimeSeries.knit 1 サンプルデータ作成 2 日付単位に集計する 2.1 月毎集計 2.2 四半期毎集計 ...