File:MagaplotFragment.jpg

#include<math.h>
#include<stdio.h>
#include <stdlib.h>
#include <complex>
//using namespace std;
#define z_type std::complex<double>
#define Re(x) x.real()
#define Im(x) x.imag()
#define RI(x) x.real(),x.imag()
#define DB double
#define DO(x,y) for(x=0;x<y;x++)

void fft(z_type *a, int N, int o) // Arry is FIRST!
{z_type u,w,t; int i,j,k,l,e=1,L,p,q,m;
q=N/2; p=2; for(m=1;p<N;m++) p*=2;
p=N-1; z_type y=z_type(0.,M_PI*o); j=0;
for(i=0;i<p;i++){ if(i<j) { t=a[j]; a[j]=a[i]; a[i]=t;}
                   k=q; while(k<=j){ j-=k; k/=2;}
                   j+=k; }
for(l=0;l<m;l++)
{ L=e; e*=2; u=1.; w=exp(y/((double)L));
  for(j=0;j<L;j++)
  { for(i=j;i<N;i+=e){k=i+L; t=a[k]*u; a[k]=a[i]-t; a[i]+=t;}
     u*=w;
} }
}

void ado(FILE *O, int X, int Y)
{ fprintf(O,"%c!PS-Adobe-2.0 EPSF-2.0\n",'%');
        fprintf(O,"%c%cBoundingBox: 0 0 %d %d\n",'%','%',X,Y);
        fprintf(O,"/M {moveto} bind def\n");
        fprintf(O,"/L {lineto} bind def\n");
        fprintf(O,"/S {stroke} bind def\n");
        fprintf(O,"/s {show newpath} bind def\n");
        fprintf(O,"/C {closepath} bind def\n");
        fprintf(O,"/F {fill} bind def\n");
        fprintf(O,"/o {.01 0 360 arc C S} bind def\n");
        fprintf(O,"/times-Roman findfont 20 scalefont setfont\n");
        fprintf(O,"/W {setlinewidth} bind def\n");
        fprintf(O,"/RGB {setrgbcolor} bind def\n");}

// #include "ado.cin"
// #include"fafo.cin"
//#include "mori.cin"
DB L1= 2.404825557695773;
DB L2= 5.5200781102863115;
DB L3= 8.653727912911013;
DB morin(DB x){ return j0(L1*x)/(1-x*x);}
DB mori0(DB x){ int n,m; DB s, xx=x*x;
DB c[16]={ 1., -0.4457964907366961303, 0.07678538241994023453, -0.0071642885058902232688,
0.00042159522055140947688, -0.000017110542281627483109, 5.0832583976057607495e-7, -1.1537378620148452816e-8,
2.0662789231930073316e-10, -2.9948657413756059965e-12, 3.5852738451127332173e-14,-3.6050239634659700777e-16,
3.0877184831292878827e-18, -2.2798156440952688462e-20, 1.4660907878585489441e-22,-8.2852774398657968065e-25};
// 16th term seems to fail; perhaps, due to the C++ rounding errors.
//with m=15, at |x|<2, the error is of order of 10^(-16)
//In this sense, the result is accurate while |x|<2.
m=15; s=c[m]*xx; for(n=m-1;n>0;n--){ s+=c[n]; s*=xx;}
return 1.+s;}
DB mori(DB x){if(fabs(x)<2) return mori0(x);
                        return morin(x);}

DB fit(DB x){ DB a,b,c,cc,s; a=-.055; b=0.02; c=0.45; cc=c*c;
s=x*(a+x*(b+x*cc)); return c*x*sqrt(x/(1.+s)); }

DB fit2(DB x){ DB q; q=5.+x*(-.5+x*(.25+x)); return x*sqrt(x/q);}

int main(){z_type * a, *b, c; int j,m,n,N; FILE *o; DB scale,step,x,y,q,p,u,v;
n=15;
N=pow(2,n);
//scale=.5;
step=sqrt(2*M_PI/N);
scale=100*step;
DB dx=step*scale;
DB dp=step/scale;
printf("2^%2d=%8d scale=%6.4lf dx=%9.8lf dp=%9.8lf\n",n,N,scale,dx,dp);
a=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));
b=(z_type *) malloc((size_t)((N+1)*sizeof(z_type)));

DO(m,N){x=m*dx; y=mori(sqrt(x)); y*=y; a[m]=y; b[m]=y*dx; }
b[0]*=.5;
fft(b,N,1); //DO(j,N) printf("%2d %18.15f %18.15f %18.15f %18.15f\n", j, RI(a[j]), RI(b[j]) );
//o=fopen("32.eps","w"); ado(o,1008,228);
o=fopen("magaplo.eps","w"); ado(o,1008,208);
#define M(x,y) fprintf(o,"%6.4f %6.4f M\n",0.+x,0.+y);
#define L(x,y) fprintf(o,"%6.4f %6.4f L\n",0.+x,0.+y);
#define o(x,y) fprintf(o,"%6.4f %6.4f o\n",0.+x,0.+y);
fprintf(o,"4 104 translate 100 100 scale\n 1 setlinejoin 2 setlinecap\n");
for(m=0;m<11;m++) {M(m,-1) L(m,1)}
for(n=-1;n<2;n++) {M(0,n) L(10,n)} fprintf(o,".002 W S\n");

DO(m,N){x=dx*m; y=Re(a[m]); if(m==0)M(x,y)else L(x,y);if(x>10) break;} fprintf(o,".008 W 0 .8 0 RGB S\n");
//DO(n,N){p=dp*n; y=fit(p); if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".008 W 1 0 0 RGB S\n");
//DO(n,N){p=dp*n; y=fit2(p); if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".008 W 0 0 1 RGB S\n");
fprintf(o,".01 W 0 0 1 RGB\n");
for(n=4;n<N;n+=4){p=dp*n; y=fit2(p); o(p,y);if(p>10) break;}

DO(n,N){p=dp*n; y=Re(b[n]); if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".02 W 0 0 1 RGB S\n");
DO(n,N){p=dp*n; y=Im(b[n]); if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".02 W 1 0 0 RGB S\n");
DO(n,N){p=dp*n; u=Re(b[n]); v=Im(b[n]); y=1-u*u-v*v; if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".02 W 0 0 0 RGB S\n");
// fprintf(o,".01 W 0 0 0 RGB\n");
//for(n=4;n<N;n+=4){p=dp*n; u=Re(b[n]); v=Im(b[n]); y=1-u*u-v*v; o(p,y);if(p>10) break;}

//DO(n,N){p=dp*n; u=Re(b[n]); v=Im(b[n]); y=fit(p)-(1.-u*u-v*v); y*=50; if(n==0)M(p,y)else L(p,y);if(p>10) break;} fprintf(o,".004 W .8 0 0 RGB S\n");
DO(n,N){p=dp*n; u=Re(b[n]); v=Im(b[n]); y=fit2(p)-(1.-u*u-v*v); y*=50; if(n==0)M(p,y)else L(p,y);if(p>10) break;}
fprintf(o,".008 W 0 0 1 RGB S\n");

fprintf(o,"showpage\n%cTrailer",'%'); fclose(o);
 system("epstopdf magaplo.eps");
 system( "open magaplo.pdf");
 free(a);
 free(b);
}

\documentclass[12pt]{article}
\usepackage{geometry}
\usepackage{graphics}
\usepackage{rotating}
\paperwidth 420pt
\paperheight 216pt
\textwidth 1420pt
\textheight 300pt
\topmargin -108pt
\oddsidemargin -73pt
\newcommand \ds {\displaystyle}
\newcommand \sx {\scalebox}
\newcommand \rme {\mathrm{e}}
\newcommand \rot {\begin{rotate}}
\newcommand \ero {\end{rotate}}
\newcommand \ing {\includegraphics}
\pagestyle{empty}
\parindent 0pt
\begin{document}
\begin{picture}(410,214)
\put(10,0){\ing{magaplo}}
%\put(.6,120){\sx{1.25}{$y$}}
\put(2,199){\sx{1.3}{$1$}}
%\put(.6,56){\sx{1.4}{$\frac{1}{2}$}}
\put(2, 100){\sx{1.3}{$0$}}
\put(-2, 0){\sx{1.3}{$-\!1$}}
\put(111,88){\sx{1.33}{$1$}}
\put(211,88){\sx{1.33}{$2$}}
\put(311,89){\sx{1.33}{$3$}}
\put(406,89){\sx{1.4}{$x$}}
%\put(411,89){\sx{1.33}{$4$}}
%\put(512,89){\sx{1.33}{$5$}}
%\put(613,89){\sx{1.33}{$6$}}
%\put(713,89){\sx{1.33}{$7$}}
%\put(813,89){\sx{1.33}{$8$}}
%\put(914,89){\sx{1.33}{$9$}}
%\put(1010,89){\sx{1.4}{$x$}}
%\put(15,192){\rot{-36}\sx{1.3}{$y\!=\! \mathrm{nori}(x)$}\ero}

\put(280,199){\rot{3}\sx{1.3}{$y\!=\! \mathrm{fit2}(x))$}\ero}
\put(280,178){\rot{3}\sx{1.3}{$y\!=\! \mathrm{maga}(x))$}\ero}

\put(118,130){\rot{-16}\sx{1.3}{$y\!=\! \mathrm{nori}(x)$}\ero}
\put(292,142){\rot{-4}\sx{1.3}{$y\!=\! \Im(\mathrm{naga}(x))$}\ero}
\put(292,120){\rot{-4}\sx{1.3}{$y\!=\! \Re(\mathrm{naga}(x))$}\ero}

\put(170,26){\sx{1.3}{$y\!=\! 50 \Big(\mathrm{fit2}(x)-\mathrm{maga}(x)\Big)$}}
\end{picture}\end{document}