C++ Figury Chladniego

// NOTE: compile with g++ filename.cpp -std=c++11

#include <iostream>
#include <cmath>
#define DIM 512
#define DM1 (DIM-1)
#define _sq(x) ((x)*(x))                           // square
#define _cb(x) abs((x)*(x)*(x))                    // absolute value of cube
#define _cr(x) (unsigned char)(pow((x),1.0/3.0))   // cube root

double max=0;

unsigned char RD(int i,int j){
   double
//
// PARAMETRY FUNKCJI
//
      a = M_PI/4,  // Kat
      m = 3,       // M
      n = 4,       // N

      kontrast = 255, // kontrast wygenerowanego obrazu

      f=0,ca=0,sa=0,pdm=0,pdn=0,oo=0;
   ca = cos(a);
   sa = sin(a);
   pdm = M_PI*100/DIM*m;
   pdn = M_PI*100/DIM*n;
   f = sin(pdm*i*0.01)*sin(pdn*j*0.01)*ca + sin(pdn*i*0.01)*sin(pdm*j*0.01)*sa;

   oo = abs(f*kontrast);

   if(oo > max)
       max = oo;

   /* Teoretycznie `f` zwraca wartosci od -2 do +2. W praktyce wychodzi mniejszy
      zakres, najczesciej -1 - 1. Dlatego mnoze przez `f` przez `kontrast` aby poprawic kontrast obrazu
      i obcinam nadmiar. Efekt jest zadowalajacy.
   */
   if(oo > 255)
      oo = 255;

   return oo;
}

void pixel_write(int,int);
FILE *fp;

int main(){
    fp = fopen("Chladni.ppm","wb");
    fprintf(fp, "P6\n%d %d\n255\n", DIM, DIM);
    for(int j=0;j<DIM;j++)
        for(int i=0;i<DIM;i++)
            pixel_write(i,j);
    fclose(fp);
    std::cout << max << std::endl;
    return 0;
}

void pixel_write(int i, int j){
    static unsigned char color[3];
    color[0] = RD(i,j)&255;
    color[1] = color[0];
    color[2] = color[0];
    fwrite(color, 1, 3, fp);
}