File:Morimap.jpg

#include <stdio.h>
#include <stdlib.h>
#define DB double
#define DO(x,y) for(x=0;x<y;x++)
//using namespace std;
#include <complex>
typedef std::complex<double> z_type;
#define Re(x) x.real()
#define Im(x) x.imag()
#define I z_type(0.,1.)
#include "conto.cin"
#include "besselj0.cin"

DB L1= 2.404825557695773;
DB L2= 5.5200781102863115;
DB L3= 8.653727912911013;
DB L4=11.791534439014281;
DB morin(DB x){ return j0(L1*x)/(1.-x*x);} // naive representation fails at x=1.
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 relative 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);}

z_type morin(z_type x){ return BesselJ0(L1*x)/(1.-x*x);} // naive representation fails at x=1.

z_type mori0(z_type x){ int n,m; z_type 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 relative 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;}
z_type mori(z_type x){if(abs(x)<2.) return mori0(x);
                        return morin(x);}

int main(){ int j,k,m,n; DB x,y, p,q, t; z_type z,c,d;
int M=801,M1=M+1;
int N=401,N1=N+1;
DB X[M1],Y[N1], g[M1*N1],f[M1*N1], w[M1*N1]; // w is working array.
char v[M1*N1]; // v is working array
//FILE *o;o=fopen("besselj0map1.eps","w");ado(o,162,82);
//FILE *o;o=fopen("01.eps","w");ado(o,1620,820);
FILE *o;o=fopen("morima.eps","w");ado(o,1620,820);
fprintf(o,"810 410 translate\n 100 100 scale\n");
DO(m,400) X[m]=-8.+.02*m;
     X[400]=-.001;
     X[401]= .001;
for(m=402;m<M1;m++) X[m]=-8.+.02*(m-1.);
DO(n,200)Y[n]=-4.+.02*n;
      Y[200]=-.001;
      Y[201]= .001;
for(n=202;n<N1;n++) Y[n]=-4.+.02*(n-1.);
//DO(m,M1)X[m]=Y[m];
for(m=-8;m<9;m++){if(m==0){M(m,-4.1)L(m,4.1)} else{M(m,-4)L(m,4)}}
for(n=-4;n<5;n++){ M( -8,n)L(8,n)}
fprintf(o,".01 W 0 0 0 RGB S\n");
DO(m,M1)DO(n,N1){g[m*N1+n]=99999; f[m*N1+n]=99999;}
DO(m,M1){x=X[m]; //printf("%5.2f\n",x);
DO(n,N1){y=Y[n]; z=z_type(x,y);
// c=BesselJ0(z);
      c=mori(z);
     p=Re(c); q=Im(c);
     if(p>-999. && p<999. && q>-999. && q<999
     )
{g[m*N1+n]=p;
 f[m*N1+n]=q;
}
                     }}
//#include "plodi.cin"
fprintf(o,"1 setlinejoin 1 setlinecap\n"); p=4;q=.2;
for(m=-4;m<4;m++)for(n=1;n<10;n+=1)conto(o,f,w,v,X,Y,M,N,(m+.1*n),-q, q); fprintf(o,".01 W 0 .6 0 RGB S\n");
for(m=0;m<4;m++) for(n=1;n<10;n+=1)conto(o,g,w,v,X,Y,M,N,-(m+.1*n),-q, q); fprintf(o,".01 W .9 0 0 RGB S\n");
for(m=0;m<4;m++) for(n=1;n<10;n+=1)conto(o,g,w,v,X,Y,M,N, (m+.1*n),-q, q); fprintf(o,".01 W 0 0 .9 RGB S\n");
for(m=1;m<9;m++) conto(o,f,w,v,X,Y,M,N, (0.-m),-p,p); fprintf(o,".02 W .9 0 0 RGB S\n");
for(m=1;m<9;m++) conto(o,f,w,v,X,Y,M,N, (0.+m),-p,p); fprintf(o,".02 W 0 0 .9 RGB S\n");
                conto(o,f,w,v,X,Y,M,N, (0. ),-2*p,2*p); fprintf(o,".02 W .6 0 .6 RGB S\n");
for(m=-8;m<0;m++) conto(o,g,w,v,X,Y,M,N, (0.+m),-p,p); fprintf(o,".02 W 0 0 0 RGB S\n");
    m=0; conto(o,g,w,v,X,Y,M,N, (0.+m),-2*p,2*p); fprintf(o,".02 W 0 0 0 RGB S\n");
for(m=1;m<9;m++) conto(o,g,w,v,X,Y,M,N, (0.+m),-p,p); fprintf(o,".02 W 0 0 0 RGB S\n");
//#include "plofu.cin"
fprintf(o,"showpage\n%c%cTrailer",'%','%'); fclose(o);
     system("epstopdf morima.eps");
     system( "open morima.pdf");
}

\documentclass[12pt]{article}
\paperheight 838px
\paperwidth 1644px
\textwidth 1294px
\textheight 1200px
\topmargin -80px
\oddsidemargin -80px
\usepackage{graphics}
\usepackage{rotating}
\usepackage{color}%<br>
\newcommand \sx {\scalebox}
\newcommand \rot {\begin{rotate}}
\newcommand \ero {\end{rotate}}
\newcommand \ing {\includegraphics}
\newcommand \rmi {\mathrm{i}}
\begin{document}
\newcommand \zoomax {
\put(16,820){\sx{4.4}{$y$}}
\put(16,630){\sx{4}{$2$}}
\put(16,430){\sx{4}{$0$}}
\put(-4, 230){\sx{4}{$-\!2$}}
\put(220, 5){\sx{4}{$-\!6$}}
\put(420, 5){\sx{4}{$-\!4$}}
\put(620, 5){\sx{4}{$-\!2$}}
\put(843, 5){\sx{4}{$0$}}
\put(1043, 5){\sx{4}{$2$}}
\put(1243, 5){\sx{4}{$4$}}
\put(1443, 5){\sx{4}{$6$}}
\put(1631,6){\sx{4}{$x$}}
}
\parindent 0pt
\begin{picture}(1616,816)
%\put(40,30){\sx{10}{\ing{besselj0map1}}}
%\put(40,30){\ing{02}}
\put(40,30){\ing{morima}}
\zoomax
\put(170,430){\sx{4}{\rot{0}$v\!=\!0$\ero}}
\put(865,470){\sx{4}{\rot{90}$v\!=\!0$\ero}}
\put(872,442){\sx{4}{\rot{50}$u\!=\!1$\ero}}
\put(862,420){\sx{4}{\rot{-51}$u\!=\!1$\ero}}
\put(700,445){\sx{3.6}{\rot{0}$v\!=\!0.1$\ero}}
%\put(660,430){\sx{4}{\rot{0}$v\!=\!0$\ero}}
\put(670,410){\sx{3.6}{\rot{0}$v\!=\!-0.1$\ero}}
\put(1092,390){\sx{4}{\rot{90}$u\!=\!0$\ero}}
\put(1144,390){\sx{4}{\rot{90}$v\!=\!0$\ero}}
\put(1224,390){\sx{4}{\rot{90}$u\!=\!0$\ero}}
\put(1280,390){\sx{4}{\rot{90}$v\!=\!0$\ero}}
\put(1354,390){\sx{4}{\rot{90}$u\!=\!0$\ero}}
\put(1414,390){\sx{4}{\rot{90}$v\!=\!0$\ero}}
\put(1488,390){\sx{4}{\rot{90}$u\!=\!0$\ero}}
\put(1548,390){\sx{4}{\rot{90}$v\!=\!0$\ero}}
%
\put(670,302){\sx{4}{\rot{8}$v\!=\!-1$\ero}}
\put(670,258){\sx{4}{\rot{8}$v\!=\!-2$\ero}}
%
\put(898,310){\sx{4}{\rot{8}$v\!=\!1$\ero}}
\put(898,262){\sx{4}{\rot{8}$v\!=\!2$\ero}}
%
\put(1026,356){\sx{3.4}{\rot{66}$u\!=\!0.1$\ero}}
%
\put(844,32){\sx{3.3}{\rot{98}$v\!=\!-8$\ero}}
\put(883,32){\sx{3.3}{\rot{80}$v\!=\!8$\ero}}
\put(929,32){\sx{3.3}{\rot{90}$u\!=\!8$\ero}}
\put(964,32){\sx{3.3}{\rot{70}$u\!=\!-8$\ero}}
\end{picture}
\end{document}