User:Kghose/ikeda code

From Wikipedia, the free encyclopedia

I used the following code on a linux box to generate the Ikeda map attractor plots and animations.

[edit] Movie Script

Use the following script to create the animation and plots. It requires ImageMagick and mencoder from mplayer.

#Notes:
#For hint about depth 8 see http://lists.mplayerhq.hu/pipermail/mencoder-users/2006-August/003817.html
#For hint about adding text using ImageMagick see http://www-128.ibm.com/developerworks/library/l-graf/?ca=dnt-428
#For bash scripting see http://linuxgazette.net/issue18/bash.html
#For mencoder options see http://www.mplayerhq.hu/DOCS/HTML/en/menc-feat-enc-images.html

#run calculations
./ikeda  

#annotate and convert to depth 8
for img in `ls ikeda0*.png`
do
  theend=${img#*0}
  caption="u = 0.${theend%.png}"
  echo $caption
  convert -font helvetica -fill black -pointsize 12 -draw "text 10,20 \"$caption\"" $img -depth 8 $img
done

#Now run mencoder on the depth 8 files
mencoder mf://ikeda0*.png -mf w=400:h=400 -ovc lavc -lavcopts vcodec=xvid -of avi -o ikeda_map.avi

[edit] Makefile

Compile the c++ code with the following makefile:

CC=g++
CFLAGS= -c -Wall -O4 -DNO_FREETYPE
LDFLAGS= -lfreetype -lm -lpng -lpngwriter -lz -DNO_FREETYPE -O4
EXECUTABLE = ikeda 

all: $(HEADERS) $(SOURCES) $(EXECUTABLE)
        
ikeda: main.cpp
        $(CC) $(LDFLAGS) main.cpp pngwriter.o -oikeda 

clean:
        rm -f *.o ikeda 

[edit] C++ Code

The c++ program requires pngwriter, libpng and zlib.

#include <stdlib.h>
#include <iostream>
#include <pngwriter.h>

using namespace std;

//N is the number of trajectory points
//X,Y are arrays (pre assigned) of length N that contain the trajectory
//bb_XXXX contains the bounding box of the plot
struct IkedaTrajectory
{
      int N; 
      double *X, *Y, bb_xmin, bb_xmax, bb_ymin, bb_ymax ;
};

//a convenience structure for the plot bounding box
struct Scale
{
      double bb_xmin, bb_xmax, bb_ymin, bb_ymax ;
};

//convenience structure to carry the simulation params
struct Params
{
    double u;
    int N,          // simulation steps
        N_start_points, //points in trajectory
        Nx, Ny,     //size of image
        last_N_points; //for attractor plot, how many points from a trajectory do we plot
};

struct ImageData
{
    char filename[80];
    int *X, *Y;
};

//u is the Ikeda parameter
//x,y is the initial point
//it is the ikeda trajectory structure we return our answer in
//remember to initialise it
void ikeda_trajectory(double u, double x, double y, IkedaTrajectory &it)
{
     it.X[0] = x ;
     it.Y[0] = y ;
     double t, x1 = x, y1 = y;
     for(int n = 1 ; n < it.N ; n++)
     {   
         /*
         if( it.bb_xmin > x1 ) it.bb_xmin = x1 ;
         if( it.bb_xmax < x1 ) it.bb_xmax = x1 ;
         if( it.bb_ymin > y1 ) it.bb_ymin = y1 ;
         if( it.bb_ymax < y1 ) it.bb_ymax = y1 ;
         */
                     
         t = 0.4 - 6.0/(1.0 + x1*x1 + y1*y1);
         it.X[n] = 1.0 + u*(x1*cos(t) - y1*sin(t)) ;
         it.Y[n] = u*(x1*sin(t) + y1*cos(t)) ;
         x1 = it.X[n];
         y1 = it.Y[n];
     }
}

//scale the plot data to the canvas
//Pre allocate X and Y
void scale_plot(IkedaTrajectory &it, Scale &sc, int Nx, int Ny, int *X, int *Y)
{          
     double scaleX, translateX, scaleY, translateY ;
     scaleX = (double)Nx / (sc.bb_xmax - sc.bb_xmin) ;
     translateX = sc.bb_xmin ;
     scaleY = (double)Ny / (sc.bb_ymax - sc.bb_ymin) ;
     translateY = sc.bb_ymin ;                
                                
     for(int n = 0 ; n < it.N ; n++)
     {
             X[n] = (int)( (it.X[n] - translateX) * scaleX );
             Y[n] = (int)( (it.Y[n] - translateY) * scaleY );
     }
}


void save_plot(Params p, IkedaTrajectory &it, Scale &it_scale, ImageData &imagedata)
{
    pngwriter png(p.Nx, p.Ny, 1.0, imagedata.filename);
    double r, b, x, y ; 
            
    cout << "Computing u = " << p.u << flush ;
    for(int n = 0 ; n < p.N_start_points ; n++)
    {
        r = rand()/((double)RAND_MAX + 1);
        b = rand()/((double)RAND_MAX + 1);    
        x = ( r - 0.5 ) * 20.0;
        y = ( r - 0.5 ) * 20.0;    
               
        ikeda_trajectory(p.u, x, y, it) ;
        scale_plot(it, it_scale, p.Nx, p.Ny, imagedata.X, imagedata.Y) ;
                
        for(int m = it.N - p.last_N_points ; m < it.N ; m++)
        {
                png.filledcircle_blend( imagedata.X[m], imagedata.Y[m], 2.0, 0.6, r, 0.5, b);
                png.plot( imagedata.X[m], imagedata.Y[m], r, 0.5, b );
        }
    }    
  
    cout << " --- Saving " << flush ;
    //png.setcompressionlevel(1); //messing with this does not speed up the process
    png.close();
    cout << " --- Done ! " << endl ;
}

int main(int argc, char *argv[])
{    
    Params p ;
    ImageData imagedata ;
    
    p.N = 500 ;
    p.N_start_points = 20000 ;
    p.Nx = 400 ;
    p.Ny = 400 ;
    p.last_N_points = 20 ;
                
    IkedaTrajectory it ;
    Scale it_scale ;
    
    it_scale.bb_xmax = 5;
    it_scale.bb_ymax = 4.5;
    it_scale.bb_xmin = -1;
    it_scale.bb_ymin = -1.5;
    
    
    it.N = p.N ;
    it.X = new double [it.N];
    it.Y = new double [it.N];
    imagedata.X = new int [it.N];
    imagedata.Y = new int [it.N];        
 
    for(double u = 0.0 ; u < 1.0 ; u += 0.001) 
    {
        p.u = u;
        sprintf(imagedata.filename, "ikeda%04d.png", (int)(p.u * 1000 +.5)); 
        save_plot(p, it, it_scale, imagedata);
    }
    
    delete[] it.X ;
    delete[] it.Y ;
    
    delete[] imagedata.X ;
    delete[] imagedata.Y ;
    
    return EXIT_SUCCESS;
}