bd6d33fb0c406a87254dfec61434e946b9c91e3f
[scilab.git] / scilab / modules / graphics / demos / cmplxfunc / MacCmplx.sci
1 //
2 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
3 // Copyright (C) 2001 - Bruno PINCON
4 // Copyright (C) 2005 - INRIA - Pierre MARECHAL <pierre.marechal@inria.fr>
5 //
6 // This file is distributed under the same license as the Scilab package.
7 //
8
9 function [xr,yr,zr,xi,yi,zi] = CmplxFacets(R,e,TypeDomain,TypeCut,n,StrFunc)
10
11     //  A function to compute the facets for drawing a complex function
12     //  on a square or a disk with branch(es) cut(s) on Ox or Oy
13     //
14     //  TypeDomain : "Square" or "Disk"
15     //     TypeCut : "Ox" or "Oy"
16     //           R : length of half a side of the square or radius of the disk
17     //           e : thin layer to avoid the branch(es) cut(s)
18     //           n : a scalar (for Square) or a 2-vector = [ntheta, nr]
19     //               (for Disk) for discretization
20     //     StrFunc : the string which names the complex function (this is
21     //               because primitive don't pass as function argument)
22     //
23     // Bruno (10/10/2001): macros for the complex function dem
24
25     if TypeDomain == "Square" then
26         if TypeCut == "Ox" then
27             x1 = linspace(-R, R, n);
28             y1 = linspace( e, R, floor(n/2));
29         else  // for TypeCut = "Oy" ...
30             x1 = linspace( e, R, floor(n/2));
31             y1 = linspace(-R, R, n);
32         end
33         X1 = ones(y1')*x1 ; Y1 = y1'*ones(x1);
34     else // for TypeDomain = "Disk"
35         r = linspace(0,R, n(2));
36         if TypeCut == "Ox" then
37             theta = linspace(0,%pi,n(1))';
38             X1 = cos(theta)*r;
39             Y1 = e + sin(theta)*r;
40         else // for TypeCut = "Oy"
41             theta = linspace(-%pi/2,%pi/2,n(1))';
42             X1 = e + cos(theta)*r;
43             Y1 = sin(theta)*r;
44         end
45     end
46     X2 = -X1 ; Y2 = -Y1;
47     Z1 = evstr(StrFunc+"(X1 + %i*Y1)");
48     Z2 = evstr(StrFunc+"(X2 + %i*Y2)");
49     [xr1,yr1,zr1] = nf3d(X1,Y1,real(Z1));
50     [xr2,yr2,zr2] = nf3d(X2,Y2,real(Z2));
51     xr = [xr1 xr2]; yr = [yr1 yr2]; zr = [zr1 zr2];
52     [xi1,yi1,zi1] = nf3d(X1,Y1,imag(Z1));
53     [xi2,yi2,zi2] = nf3d(X2,Y2,imag(Z2));
54     xi = [xi1 xi2]; yi = [yi1 yi2]; zi = [zi1 zi2];
55 endfunction
56
57
58 function []=PlotCmplxFunc(R,e,TypeDomain,TypeCut,n,StrFunc,theta,alpha,DomReal)
59
60     //  A function to draw on a square or a disk a complex function
61     //  with branch(es) cut(s) on Ox or Oy
62     //
63     //  TypeDomain : "Square" or "Disk"
64     //     TypeCut : "Ox" or "Oy"
65     //           R : length of half a side of the square or radius of the disk
66     //           e : thin layer to avoid the branch(es) cut(s)
67     //           n : a scalar (for Square) or a 2-vector = [ntheta, nr]
68     //               (for Disk) for discretization
69     //     StrFunc : the string which names the complex function (this is
70     //               because primitive don't pass as function argument)
71     // theta,alpha : usual parameters for plot3d
72     //     DomReal : interval for which the real restriction is drawn
73
74     // Bruno (10/10/2001): macros for the complex function dem
75
76     // Adapted for new graphic by Pierre MARECHAL ( 16/12/2005 )
77
78     // computes the facets
79
80     [xr,yr,zr,xi,yi,zi] = CmplxFacets(R,e,TypeDomain,TypeCut,n,StrFunc)
81
82     // draw
83     // ============================================
84
85     current_figure_hdl = scf(100001);
86     clf(current_figure_hdl,"reset");
87
88     drawlater();
89
90     // Title
91     // ============================================
92
93     my_title_axes = newaxes();
94
95     // make axes transparent
96     my_title_axes.filled = "off";
97
98     Rs = string(R);
99
100     if TypeDomain == "Square" then
101         end_title = " Function on [-"+Rs+","+Rs+"]x[-"+Rs+","+Rs+"]"
102     else
103         end_title = " Function on D(0,R="+Rs+")"
104     end
105
106     if StrFunc == "f" then
107         the_title = "Your Custom (named f) Complex" + end_title;
108     else
109         the_title = "The Complex " + StrFunc + end_title;
110     end
111
112     xtitle(the_title);
113
114     my_title_axes.title.text       = the_title;
115     my_title_axes.title.font_size  = 3;
116     my_title_axes.title.font_style = 2;
117     my_title_axes.margins = [ 0.08 0.08 0.08 0.08 ]
118
119     // plot Im(z)
120     // ============================================
121
122     subplot(1,2,1);
123     plot3d(xi,yi,zi,theta,alpha,"Re(z)@Im(z)@",[2 6 4]);
124     xtitle("Im("+StrFunc+"(z))");
125
126     // plot Re(z) + the real restriction
127     // ============================================
128
129     subplot(1,2,2);
130     plot3d(xr,yr,zr,theta,alpha,"Re(z)@Im(z)@",[ 2 6 4]);
131     xtitle("Re("+StrFunc+"(z))");
132
133     // real function in yellow
134     // ============================================
135
136     if DomReal(2) > DomReal(1) then
137         //xstring(0.1,-0.15," In yellow : the real "+StrFunc+" function")
138         xx = linspace(DomReal(1),DomReal(2),40)';
139         yy = zeros(xx);
140         zz = evstr(StrFunc+"(xx)");
141         param3d1(xx,yy,list(zz,32),theta,alpha,flag=[0,0]);
142         yellow_line = get("hdl");
143         yellow_line.thickness = 3;
144         captions(yellow_line, "the real "+StrFunc+" function", "lower_caption");
145     end
146
147
148
149     drawnow();
150
151 endfunction