Processingで複素関数を可視化する

Processingとはアートやデザインなど視覚的表現を簡単に行えるプログラミング言語です。
今回はこのProcessingを使って複素関数の可視化を行ってみたいと思います。

## この記事のゴール
以下のように、三角関数sin(z)を複素平面上で可視化します。
![output_sinz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010852/c95a7308-ef4b-486d-2ef5-bdc56d4631c1.png)


## 基本方針
複素関数は２次元の入力に対して2次元の出力が行われるため、可視化しようとすると4次元空間が必要になってしまいます。
そこで「定義域の着色」という手法を用いて可視化を行います。
（詳しくは[こちらの記事](https://sw1227.hatenablog.com/entry/2018/12/04/090000#13-Domain-Coloring%E3%81%AB%E3%82%88%E3%82%8B%E8%A7%A3%E6%B1%BA)がわかりやかったです）

## さっそくやってみる
今回は三角関数sin(z)にz=x+iyを入力したときの様子を可視化してみたいと思います。
タイトルに複素関数を可視化すると書きましたが、今回は汎用的な方法ではなくあくまでsin(z)の可視化を行う方法のみを解説します。

## 事前準備
Processingのセットアップを行います。

```java
void setup()
{
  size(1000,1000);
  noStroke();
  colorMode(HSB, 1, 1, 1);
  loadPixels();
  drawComplexFunction();
  updatePixels();
  save("output.png");
}
```
drawComplexFunction();に実際の処理を書いていきます。

## 複素数のクラスを用意
複素数を扱うために、簡単な複素数のクラスを用意します。

```java
class Complex
{
  public float x,y;
  
  Complex(float x_, float y_)
  {
    x = x_;
    y = y_;
  }
  
  float absolute()
  {
    return sqrt(x*x + y*y);
  }
  
  float arg()
  {
    if(x>0){
      return atan(y/x);
    }
    if(x<0 && y>=0){
      return atan(y/x) + PI;
    }
    if(x<0 && y<0){
      return atan(y/x) - PI;
    }
    if(x==0 && y>0){
      return HALF_PI;
    }
    if(x==0 && y<0){
      return -HALF_PI;
    }
    else{
      return 0;
    }      
  }
}
```
複素数の偏角は次式で計算できるようなので、この通りに実装しました。
![スクリーンショット 2021-09-12 11.41.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010852/647a4249-b913-f3d8-57ba-a7a44a0c3206.png)


## 各画素ごとに、色を塗る
Processingのキャンバスを複素平面に見立て、各ピクセル(x,y)をz=x+iyに対応させ、zを三角関数sin(z)に入力したときの出力z'を算出します。
そしてz'の偏角arg(z')を色相、絶対値|z'|を彩度に割り当てて可視化します。
この時点ではまたsin(z)は実装していません。

```java
void drawComplexFunction()
{
  for(int x = -width/2; x < width/2; x += 1)
  {
    for(int y = -height/2; y < height/2; y += 1)
    {
      float x_ = x * PI / (width/2);
      float y_ = y * PI / (height/2);
      Complex z = new Complex(x_, y_);
      
      Complex z_ = sinz(z);
      
      float hue = (z_.arg() + PI) / TWO_PI;
      float saturation = z_.absolute();
      pixels[(x+width/2) + (y+height/2)*width] = color(hue, saturation, 1); 
    }
  }
}
```

## sin(z)を実装する
組み込みのsin関数では、複素数の計算ができません。
そのため、自作の関数sinzを実装します。

三角関数のマクローリン展開を利用します。
![スクリーンショット 2021-09-12 11.52.40.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010852/764c01a9-ce78-8131-0e31-026dcf611b74.png)


これをある程度の項数まで実装して、近似的な関数とします。

普通にz=x+iyを代入して計算するととてつもなく大変なので、
![スクリーンショット 2021-09-12 11.58.35.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010852/9c68d390-f518-59f2-1a5c-ddc828bec55d.png)
を利用することで次のように簡単になります（5乗、7乗…も同様）。

```java
Complex sinz(Complex z)
{
  float u = z.x
    - pow(z.absolute(), 3) * cos(3*z.arg()) / (3*2)
    + pow(z.absolute(), 5) * cos(5*z.arg()) / (5*4*3*2)
    - pow(z.absolute(), 7) * cos(7*z.arg()) / (7*6*5*4*3*2)
    + pow(z.absolute(), 9) * cos(9*z.arg()) / (9*8*7*6*5*4*3*2);
  float v = z.y
    - pow(z.absolute(), 3) * sin(3*z.arg()) / (3*2)
    + pow(z.absolute(), 5) * sin(5*z.arg()) / (5*4*3*2)
    - pow(z.absolute(), 7) * sin(7*z.arg()) / (7*6*5*4*3*2)
    + pow(z.absolute(), 9) * sin(9*z.arg()) / (9*8*7*6*5*4*3*2);
  return new Complex(u, v);
}
```

## 完成

以上でコードは完成となります。
これを実行すると以下の画像が得られます。
![output_sinz.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2010852/c95a7308-ef4b-486d-2ef5-bdc56d4631c1.png)

最終型のコードは以下となります。

```java
void setup()
{
  size(1000,1000);
  noStroke();
  colorMode(HSB, 1, 1, 1);
  loadPixels();
  drawComplexFunction();
  updatePixels();
  save("output.png");
}

void drawComplexFunction()
{
  for(int x = -width/2; x < width/2; x += 1)
  {
    for(int y = -height/2; y < height/2; y += 1)
    {
      float x_ = x * PI / (width/2);
      float y_ = y * PI / (height/2);
      Complex z = new Complex(x_, y_);
      
      Complex z_ = sinz(z);
      
      float hue = (z_.arg() + PI) / TWO_PI;
      float saturation = z_.absolute();
      pixels[(x+width/2) + (y+height/2)*width] = color(hue, saturation, 1); 
    }
  }
}

Complex sinz(Complex z)
{
  float u = z.x
    - pow(z.absolute(), 3) * cos(3*z.arg()) / (3*2)
    + pow(z.absolute(), 5) * cos(5*z.arg()) / (5*4*3*2)
    - pow(z.absolute(), 7) * cos(7*z.arg()) / (7*6*5*4*3*2)
    + pow(z.absolute(), 9) * cos(9*z.arg()) / (9*8*7*6*5*4*3*2);
  float v = z.y
    - pow(z.absolute(), 3) * sin(3*z.arg()) / (3*2)
    + pow(z.absolute(), 5) * sin(5*z.arg()) / (5*4*3*2)
    - pow(z.absolute(), 7) * sin(7*z.arg()) / (7*6*5*4*3*2)
    + pow(z.absolute(), 9) * sin(9*z.arg()) / (9*8*7*6*5*4*3*2);
  return new Complex(u, v);
}

class Complex
{
  public float x,y;
  
  Complex(float x_, float y_)
  {
    x = x_;
    y = y_;
  }
  
  float absolute()
  {
    return sqrt(x*x + y*y);
  }
  
  float arg()
  {
    if(x>0){
      return atan(y/x);
    }
    if(x<0 && y>=0){
      return atan(y/x) + PI;
    }
    if(x<0 && y<0){
      return atan(y/x) - PI;
    }
    if(x==0 && y>0){
      return HALF_PI;
    }
    if(x==0 && y<0){
      return -HALF_PI;
    }
    else{
      return 0;
    }      
  }
}
```

最後までお読みいただきありがとうございました。