#include <stdio.h> 
#include <stdlib.h>
#include <string.h>
#include <math.h> 
#define Nx 500    
#define Ny 500    
float bessel1(float x);
float bessj1(float x);
/*
    single source airy disk:
    I_x=I_0 {2*J_1(x)/x}^2
    x= pi*D/lambda *sin(theta)
*/
int main() {
    float x,y,x0,y0,r,j1,Ix;
    float x1,y1;
    float rscale=.04;
    int i,j,result,array[Nx][Ny];
    int irows=Nx,icols=Ny, iframe=0;
    char namebuf[13];
    x0=150.;     //pixels=coordinates
    y0=250.;
    x1=350.;     //pixels=coordinates
    y1=250.;
    for(iframe=1;iframe<=200;iframe++){
        x0=50.+(float)iframe;
        x1=450-(float)iframe;
        for(i=0 ; i<irows ; i++) {
            for(j=0 ; j<icols ; j++) {
                x=j-x0;
                y=i-y0;
                r=pow(x*x+y*y,.5)*rscale;
                if(r<.000000001){r=.000000001;}    // watch that over/underflow
                j1=bessj1(r);
                Ix=2.*j1*j1/(r*r);
                x=j-x1;
                y=i-y1;
                r=pow(x*x+y*y,.5)*rscale;
                if(r<.000000001){r=.000000001;}    // watch that over/underflow
                j1=bessj1(r);
                Ix=Ix+2.*j1*j1/(r*r);
                array[i][j]=(int)(64000.*Ix);
                //printf("for x= %f for y= %f j1=%f Ix=%f I=%i \n",x,y,j1,Ix,array[i][j]);
            }
        }
        sprintf( namebuf, "b_hf%.4i.tga", iframe );
        result=tgaout(irows, icols, namebuf, &array[0][0]);
        printf("Writing  frame %i\n",iframe);
        sprintf( namebuf, "b_gs%.4i.tga", iframe );
        result=tgaout2(irows, icols, namebuf, &array[0][0]);
        printf("Writing  frame %i\n",iframe);
    }    

    return 0;
}    
 
//numerical recipes in C
float bessj1(float x){
    float ax,z;
    double xx,y,ans,ans1,ans2; 
    if ((ax=fabs(x)) < 8.0) { 
        y=x*x;
        ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
            +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
        ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
            +y*(99447.43394+y*(376.9991397+y*1.0))));
        ans=ans1/ans2;
    } else { 
        z=8.0/ax;
        y=z*z;
        xx=ax-2.356194491;
        ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
            +y*(0.2457520174e-5+y*(-0.240337019e-6))));
        ans2=0.04687499995+y*(-0.2002690873e-3
            +y*(0.8449199096e-5+y*(-0.88228987e-6
            +y*0.105787412e-6)));
        ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
        if (x < 0.0) ans = -ans;
    }
    return ans;
}


/* Function tgaout2
 * to write integer data as uncompressed targa (grayscale)
 * 
 * adapted by mrg3@psu.edu from
 * 
 * asc2tga.cpp
 * converts ascii matrix to tga heightfield
 * Copyright (C) 1999 by Ortwin Glueck
 * Released under the GNU General Public License (GPL)
 * This is free software. No warranty.
 */

int tgaout2(int irows, int icols, char *fname, int *data)
{
unsigned char tgaheader[]={0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,24,32}; 
 FILE *out;
 unsigned int value,i,j;
 unsigned short rows,cols;
 unsigned char BGR[3];
 unsigned short *target;
 float fvalue;

 out=fopen(fname,"wb");
 if (out==NULL)
 {
  printf("Could not open output file\n");
  exit(1);
 }
 
 target=(unsigned short*)&tgaheader[12];
 (unsigned short)rows=irows;
 (unsigned short)cols=icols;
 target[0]=rows;
 target[1]=cols;
 fwrite(tgaheader,18,1,out);
 
 BGR[0]=0; //blue is not used
 for (i=0;i<rows;i++)
  for (j=0;j<cols;j++)
  {
   value = *(data+i*icols+j);
   if(value<0){
        value=0;
   }
   if(value>65535){
        value=65535;
    }
    fvalue=(float)(value/65535.*253.+2.); //all rgb channels get this much, with a little base
   //printf("Writing %u value\n",value);
   //*((unsigned short*)&BGR[1])=value;
   BGR[0]=(unsigned char)fvalue;
   BGR[1]=BGR[0];
   BGR[2]=BGR[0];
   
   fwrite(BGR,3,1,out);
  }
    fclose(out);
    return 0;

}



/* Function tgaout
 * to write integer data to uncompressed targa heightfield
 * 
 * adapted by mrg3@psu.edu from
 * 
 * asc2tga.cpp
 * converts ascii matrix to tga heightfield
 * Copyright (C) 1999 by Ortwin Glueck
 * Released under the GNU General Public License (GPL)
 * This is free software. No warranty.
 */

int tgaout(int irows, int icols, char *fname, int *data)
{
unsigned char tgaheader[]={0,0,2,0,0,0,0,0,0,0,0,0,0,0,0,0,24,32}; 
 FILE *out;
 unsigned int value,i,j;
 unsigned short rows,cols;
 unsigned char BGR[3];
 unsigned short *target;

 out=fopen(fname,"wb");
 if (out==NULL)
 {
  printf("Could not open output file\n");
  exit(1);
 }
 
 target=(unsigned short*)&tgaheader[12];
 (unsigned short)rows=irows;
 (unsigned short)cols=icols;
 target[0]=rows;
 target[1]=cols;
 fwrite(tgaheader,18,1,out);
 
 BGR[0]=0; //blue is not used
 for (i=0;i<rows;i++)
  for (j=0;j<cols;j++)
  {
   value = *(data+i*icols+j);
   if(value<0){
        value=0;
   }
   if(value>65535){
        value=65535;
    }

   //printf("Writing %u value\n",value);
   *((unsigned short*)&BGR[1])=value;
   fwrite(BGR,3,1,out);
  }
    fclose(out);
    return 0;

}