1 /***************************************************************************
2 * cont.cpp is part of Math Graphic Library
3 * Copyright (C) 2007-2016 Alexey Balakin <mathgl.abalakin@gmail.ru> *
4 * *
5 * This program is free software; you can redistribute it and/or modify *
6 * it under the terms of the GNU Lesser General Public License as *
7 * published by the Free Software Foundation; either version 3 of the *
8 * License, or (at your option) any later version. *
9 * *
10 * This program is distributed in the hope that it will be useful, *
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13 * GNU General Public License for more details. *
14 * *
15 * You should have received a copy of the GNU Lesser General Public *
16 * License along with this program; if not, write to the *
17 * Free Software Foundation, Inc., *
18 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19 ***************************************************************************/
20 // NOTE: Borland before 2007 (i.e. < 0x0600) use <algorithm.h>, after 0x0630 it use <algorithm>.
21 // I don't find information about 2009, 2010 versions (i.e. 0x0610 and 0x0620).
22 // May be condition below can be rewritten as (__CODEGEARC__ >= 0x0600)
23 #if !defined(__BORLANDC__) || (__CODEGEARC__ >= 0x0630)
24 #include <algorithm>
25 #else
26 #include <algorithm.h>
27 #endif
28 #include <list>
29
30 #include "mgl2/surf.h"
31 #include "mgl2/cont.h"
32 #include "mgl2/data.h"
33 #include "mgl2/eval.h"
34 #include "mgl2/font.h"
35 #include "mgl2/base.h"
36 #include "cont.hpp"
37 //-----------------------------------------------------------------------------
38 //
39 // Text printing along a curve
40 //
41 //-----------------------------------------------------------------------------
mgl_string_curve(mglBase * gr,long f,long,const long * ff,const long * nn,const wchar_t * text,const char * font,mreal size)42 void MGL_NO_EXPORT mgl_string_curve(mglBase *gr,long f,long ,const long *ff,const long *nn,const wchar_t *text, const char *font, mreal size)
43 {
44 if(f<0 || nn[f]<0) return; // do nothing since there is no curve
45 if(!font) font="";
46 int pos = strchr(font,'T') ? 1:-1, align;
47 bool cc=mglGetStyle(font,0,&align); align = align&3;
48 mreal h=gr->TextHeight(font,size)/2, g = 1.1*h;
49 wchar_t L[2]=L"a";
50
51 if(align==1) // TODO divide curve by 2
52 {}
53
54 std::vector<mglPoint> qa, qb; // curves above and below original
55 if(ff[f]<0) for(long i=nn[f];i>=0 && i!=f;i=nn[i]) // find first real point
56 if(ff[i]>=0) { f=i; break; }
57 if(ff[f]<0) return;
58 mreal c=cc?gr->AddTexture(font) : gr->GetClrC(ff[f]);
59 mglPoint p=gr->GetPntP(ff[f]), q=p, s;
60 for(long i=nn[f];i>=0 && i!=f;i=nn[i]) // find second real point
61 if(ff[i]>=0) { s=gr->GetPntP(ff[i]); break; }
62 mglPoint l=!(s-q), t=l;
63 qa.push_back(q+l*g); qb.push_back(q-l*h);
64 for(long i=nn[f];i>=0 && i!=f;i=nn[i]) // construct curves
65 {
66 p=q; q=s; l=t;
67 if(nn[i]>=0 && ff[nn[i]]>=0) { s=gr->GetPntP(ff[nn[i]]); t=!(s-q); }
68 mreal tet = t.x*l.y-t.y*l.x;
69 mreal tt = 1+fabs(t.x*l.x+t.y*l.y);
70 if(tet>0)
71 { qa.push_back(q+l*g); qa.push_back(q+t*g); qb.push_back(q-(l+t)*(h/tt)); }
72 else if(tet<0)
73 { qb.push_back(q-l*h); qb.push_back(q-t*h); qa.push_back(q+(l+t)*(g/tt)); }
74 else
75 { qa.push_back(q+l*g); qb.push_back(q-l*h); }
76 }
77 if(pos>0) qa=qb;
78 // adjust text direction
79 bool rev = align==2;
80 const char *ffont=mglchr(font,':');
81 char *fnt = new char[strlen(font)+5];
82 if(ffont) strcpy(fnt,ffont); else *fnt=0;
83 /* if(qa[0].x>qa[1].x)
84 {
85 if(align==0){ strcat(fnt,":R"); align=2; }
86 else if(align==1) rev = true;
87 else { strcat(fnt,":L"); align=0; }
88 }*/
89 if(mglchr(font,'T')) strcat(fnt,":T");
90 if(rev) reverse(qa.begin(),qa.end());
91 long len = mgl_wcslen(text);
92 mreal *wdt=new mreal[len+1];
93 for(long j=0;j<len;j++) { L[0]=text[j]; wdt[j]=1.2*gr->TextWidth(L,font,size); }
94 wdt[len]=0;
95
96 // place glyphs points
97 mglPoint *pt=new mglPoint[len+1];
98 pt[0] = qa[0];
99 long i=0, k, m = qa.size();
100 mreal t1,t2, tt=0;
101 for(long j=0;j<len;j++)
102 {
103 mreal w = align==1 ? wdt[j] : (wdt[j]+wdt[j+1])/2;
104 p = pt[j]; w *= 1.03;
105 for(k=i+1;k<m;k++) if((p-qa[k]).norm()>w) break;
106 if(k>i+1 && k<m) tt=-1;
107 i = k<m ? k-1 : m-2; // check if end of curve
108 q = qa[i]; s = qa[i+1]; // points of line segment
109 mreal a = (q-s)*(q-s), b = (q-p)*(q-s), d = (q-p)*(q-p)-w*w;
110 w = sqrt(b*b-a*d); // NOTE: b*b>a*d should be here!
111 if(b*b>1e3*a*d) { t1 = d/(b+w); t2 = d/(b-w); } // keep precision
112 else { t1 = (b-w)/a; t2 = (b+w)/a; }
113 if(t1<0 || t1<tt) t1=t2; // t1<t2 should be here!
114 tt=t1; pt[j+1] = q+(s-q)*tt;
115 }
116 if(rev) pos=-pos;
117 mreal dc = (cc && len>1)?1/MGL_FEPSILON/(len-1):0;
118 for(long j=0;j<len;j++) // draw text
119 {
120 L[0] = text[align!=2?j:len-1-j]; s = pt[j+1]-pt[j]; l = !s;
121 gr->text_plot(gr->AddPnt(pt[j]+(pos*h)*l,c+dc*i,align!=2?s:-s,-1,-1),L,fnt,size,0.05,c+dc*j);
122 }
123 delete []wdt; delete []pt; delete []fnt;
124 }
125 //-----------------------------------------------------------------------------
mgl_textw_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,const wchar_t * text,const char * font,const char * opt)126 void MGL_EXPORT mgl_textw_xyz(HMGL gr, HCDT x, HCDT y, HCDT z,const wchar_t *text, const char *font, const char *opt)
127 {
128 long n=y->GetNx();
129 if(mgl_check_dim1(gr,x,y,z,0,"Text")) return;
130
131 gr->SaveState(opt);
132 static int cgid=1; gr->StartGroup("TextC",cgid++);
133
134 long kq = gr->AllocPnts(n);
135 long *nn = new long[n], *ff = new long[n];
136 #pragma omp parallel for
137 for(long i=0;i<n;i++)
138 { ff[i] = kq+i; nn[i] = i+1;
139 gr->AddPntQ(kq+i,mglPoint(x->v(i),y->v(i),z->v(i)),-1); }
140 nn[n-1]=-1;
141 mgl_string_curve(gr,0,n,ff,nn,text,font,-1);
142 delete []ff; delete []nn; gr->EndGroup();
143 }
144 //-----------------------------------------------------------------------------
mgl_textw_xy(HMGL gr,HCDT x,HCDT y,const wchar_t * text,const char * font,const char * opt)145 void MGL_EXPORT mgl_textw_xy(HMGL gr, HCDT x, HCDT y, const wchar_t *text, const char *font, const char *opt)
146 {
147 gr->SaveState(opt);
148 mglDataV z(y->GetNx()); z.Fill(gr->AdjustZMin());
149 mgl_textw_xyz(gr,x,y,&z,text,font,0);
150 }
151 //-----------------------------------------------------------------------------
mgl_textw_y(HMGL gr,HCDT y,const wchar_t * text,const char * font,const char * opt)152 void MGL_EXPORT mgl_textw_y(HMGL gr, HCDT y, const wchar_t *text, const char *font, const char *opt)
153 {
154 gr->SaveState(opt);
155 mglDataV x(y->GetNx()), z(y->GetNx());
156 x.Fill(gr->Min.x,gr->Max.x); z.Fill(gr->AdjustZMin());
157 mgl_textw_xyz(gr,&x,y,&z,text,font,0);
158 }
159 //-----------------------------------------------------------------------------
mgl_text_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,const char * text,const char * font,const char * opt)160 void MGL_EXPORT mgl_text_xyz(HMGL gr, HCDT x, HCDT y, HCDT z,const char *text, const char *font, const char *opt)
161 {
162 MGL_TO_WCS(text,mgl_textw_xyz(gr,x,y,z, wcs, font, opt));
163 }
164 //-----------------------------------------------------------------------------
mgl_text_xy(HMGL gr,HCDT x,HCDT y,const char * text,const char * font,const char * opt)165 void MGL_EXPORT mgl_text_xy(HMGL gr, HCDT x, HCDT y, const char *text, const char *font, const char *opt)
166 {
167 mglDataV z(y->GetNx()); z.Fill(gr->AdjustZMin());
168 mgl_text_xyz(gr,x,y,&z,text,font,opt);
169 }
170 //-----------------------------------------------------------------------------
mgl_text_y(HMGL gr,HCDT y,const char * text,const char * font,const char * opt)171 void MGL_EXPORT mgl_text_y(HMGL gr, HCDT y, const char *text, const char *font, const char *opt)
172 {
173 mglDataV x(y->GetNx()), z(y->GetNx());
174 x.Fill(gr->Min.x,gr->Max.x); z.Fill(gr->AdjustZMin());
175 mgl_text_xyz(gr,&x,y,&z,text,font,opt);
176 }
177 //-----------------------------------------------------------------------------
mgl_text_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,const char * text,const char * font,const char * opt,int l,int n,int lo)178 void MGL_EXPORT mgl_text_xyz_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z,const char *text,const char *font, const char *opt,int l,int n,int lo)
179 { char *s=new char[l+1]; memcpy(s,text,l); s[l]=0;
180 char *f=new char[n+1]; memcpy(f,font,n); f[n]=0;
181 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
182 mgl_text_xyz(_GR_, _DA_(x),_DA_(y), _DA_(z), s, f, o);
183 delete []o; delete []s; delete []f; }
184 //-----------------------------------------------------------------------------
mgl_text_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,const char * text,const char * font,const char * opt,int l,int n,int lo)185 void MGL_EXPORT mgl_text_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, const char *text, const char *font, const char *opt, int l,int n,int lo)
186 { char *s=new char[l+1]; memcpy(s,text,l); s[l]=0;
187 char *f=new char[n+1]; memcpy(f,font,n); f[n]=0;
188 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
189 mgl_text_xy(_GR_, _DA_(x),_DA_(y),s,f,o);
190 delete []o; delete []s; delete []f; }
191 //-----------------------------------------------------------------------------
mgl_text_y_(uintptr_t * gr,uintptr_t * y,const char * text,const char * font,const char * opt,int l,int n,int lo)192 void MGL_EXPORT mgl_text_y_(uintptr_t *gr, uintptr_t *y, const char *text, const char *font, const char *opt, int l,int n,int lo)
193 { char *s=new char[l+1]; memcpy(s,text,l); s[l]=0;
194 char *f=new char[n+1]; memcpy(f,font,n); f[n]=0;
195 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
196 mgl_text_y(_GR_, _DA_(y),s,f,o); delete []o; delete []s; delete []f; }
197 //-----------------------------------------------------------------------------
198 //
199 // DCont series
200 //
201 //-----------------------------------------------------------------------------
set(mreal u1,mreal v1,mreal u2,mreal v2,long i,long j,long k,HCDT x,HCDT y,HCDT z)202 bool mglSegment::set(mreal u1,mreal v1,mreal u2,mreal v2,long i,long j,long k,HCDT x, HCDT y, HCDT z)
203 {
204 bool res=(v1>=0 && v1<=MGL_FEPSILON && u1>=0 && u1<=MGL_FEPSILON && v2>=0 && v2<=MGL_FEPSILON && u2>=0 && u2<=MGL_FEPSILON);
205 if(v1==v2 && u1==u2) res=false; // NOTE: shouldn't be here never
206 if(res)
207 {
208 p1.Set(mgl_data_linear(x,i+u1,j+v1,k), mgl_data_linear(y,i+u1,j+v1,k), mgl_data_linear(z,i+u1,j+v1,k));
209 p2.Set(mgl_data_linear(x,i+u2,j+v2,k), mgl_data_linear(y,i+u2,j+v2,k), mgl_data_linear(z,i+u2,j+v2,k));
210 }
211 return res;
212 }
213 //-----------------------------------------------------------------------------
set(const mglPoint & q1,const mglPoint & q2,HCDT x,HCDT y,HCDT z,bool nboth)214 void mglSegment::set(const mglPoint &q1, const mglPoint &q2,HCDT x, HCDT y, HCDT z, bool nboth)
215 {
216 if(nboth)
217 {
218 p1.Set(mgl_data_linear(x,q1.x,0,0), mgl_data_linear(y,q1.y,0,0), mgl_data_linear(z,q1.z,0,0));
219 p2.Set(mgl_data_linear(x,q2.x,0,0), mgl_data_linear(y,q2.y,0,0), mgl_data_linear(z,q2.z,0,0));
220 }
221 else
222 {
223 p1.Set(mgl_data_linear(x,q1.x,q1.y,q1.z), mgl_data_linear(y,q1.x,q1.y,q1.z), mgl_data_linear(z,q1.x,q1.y,q1.z));
224 p2.Set(mgl_data_linear(x,q2.x,q2.y,q2.z), mgl_data_linear(y,q2.x,q2.y,q2.z), mgl_data_linear(z,q2.x,q2.y,q2.z));
225 }
226 }
227 //-----------------------------------------------------------------------------
mgl_dcont_add_pnt(mglPoint p1,mglPoint p2,HCDT b,mreal val,std::vector<mglPoint> & pp)228 void MGL_NO_EXPORT mgl_dcont_add_pnt(mglPoint p1, mglPoint p2, HCDT b, mreal val, std::vector<mglPoint> &pp)
229 {
230 mreal b1=b->Linear(p1.x,p1.y,p1.z), b2=b->Linear(p2.x,p2.y,p2.z);
231 mreal d = (val-b1)/(b2-b1);
232 if(d>=0 && d<=1) pp.push_back(p2*d+p1*(1-d));
233 }
234 //-----------------------------------------------------------------------------
mgl_get_dlines(mreal val,HCDT a,HCDT b,HCDT x,HCDT y,HCDT z)235 std::vector<mglSegment> MGL_EXPORT mgl_get_dlines(mreal val, HCDT a, HCDT b, HCDT x, HCDT y, HCDT z)
236 {
237 long nx=a->GetNx(), ny=a->GetNy(), nz=a->GetNz(), nn=nx*ny;
238 bool nboth = mgl_isnboth(x,y,z,a);
239 std::vector<mglSegment> lines;
240 for(long k=0;k<nz-1;k++) for(long j=0;j<ny-1;j++) for(long i=0;i<nx-1;i++)
241 {
242 std::vector<mglPoint> pp;
243 long i0 = i+nx*(j+ny*k);
244 mreal v1=a->vthr(i0), v2=a->vthr(i0+1), v3=a->vthr(i0+nx), v4=a->vthr(i0+1+nx);
245 mreal v5=a->vthr(i0+nn), v6=a->vthr(i0+1+nn), v7=a->vthr(i0+nx+nn), v8=a->vthr(i0+1+nx+nn);
246 // first find isosurface for a
247 mreal d1,d2,d3,d4,d5,d6,d7,d8,dA,dB,dC,dD;
248 d1=(val-v1)/(v2-v1); d1=(d1>=0&&d1<=1)?d1:NAN;
249 d2=(val-v2)/(v4-v2); d2=(d2>=0&&d2<=1)?d2:NAN;
250 d3=(val-v3)/(v4-v3); d3=(d3>=0&&d3<=1)?d3:NAN;
251 d4=(val-v1)/(v3-v1); d4=(d4>=0&&d4<=1)?d4:NAN;
252 mglPoint p1(i+d1,j,k), p2(i+1,j+d2,k), p3(i+d3,j+1,k), p4(i,j+d4,k);
253 d5=(val-v5)/(v6-v5); d5=(d5>=0&&d5<=1)?d5:NAN;
254 d6=(val-v6)/(v8-v6); d6=(d6>=0&&d6<=1)?d6:NAN;
255 d7=(val-v7)/(v8-v7); d7=(d7>=0&&d7<=1)?d7:NAN;
256 d8=(val-v5)/(v7-v5); d8=(d8>=0&&d8<=1)?d8:NAN;
257 mglPoint p5(i+d5,j,k+1), p6(i+1,j+d6,k+1), p7(i+d7,j+1,k+1), p8(i,j+d8,k+1);
258 dA=(val-v1)/(v5-v1); dA=(dA>=0&&dA<=1)?dA:NAN;
259 dB=(val-v2)/(v6-v2); dB=(dB>=0&&dB<=1)?dB:NAN;
260 dC=(val-v3)/(v7-v3); dC=(dC>=0&&dC<=1)?dC:NAN;
261 dD=(val-v4)/(v8-v4); dD=(dD>=0&&dD<=1)?dD:NAN;
262 mglPoint pA(i,j,k+dA), pB(i+1,j,k+dB), pC(i,j+1,k+dC), pD(i+1,j+1,k+dD);
263 // next find its cross-section with isosurface of b at faces
264 mgl_dcont_add_pnt(p1,p2,b,val,pp); mgl_dcont_add_pnt(p1,p4,b,val,pp);
265 mgl_dcont_add_pnt(p1,p3,b,val,pp); mgl_dcont_add_pnt(p1,p5,b,val,pp);
266 mgl_dcont_add_pnt(p1,pA,b,val,pp); mgl_dcont_add_pnt(p1,pB,b,val,pp);
267 mgl_dcont_add_pnt(p2,p4,b,val,pp); mgl_dcont_add_pnt(p2,p3,b,val,pp);
268 mgl_dcont_add_pnt(p2,p6,b,val,pp); mgl_dcont_add_pnt(p2,pD,b,val,pp);
269 mgl_dcont_add_pnt(p2,pB,b,val,pp); mgl_dcont_add_pnt(p4,p3,b,val,pp);
270 mgl_dcont_add_pnt(p4,p8,b,val,pp); mgl_dcont_add_pnt(p4,pA,b,val,pp);
271 mgl_dcont_add_pnt(p4,pC,b,val,pp); mgl_dcont_add_pnt(p5,p6,b,val,pp);
272 mgl_dcont_add_pnt(p5,p8,b,val,pp); mgl_dcont_add_pnt(p5,p7,b,val,pp);
273 mgl_dcont_add_pnt(p5,pA,b,val,pp); mgl_dcont_add_pnt(p5,pB,b,val,pp);
274 mgl_dcont_add_pnt(p6,p8,b,val,pp); mgl_dcont_add_pnt(p6,p7,b,val,pp);
275 mgl_dcont_add_pnt(p6,pD,b,val,pp); mgl_dcont_add_pnt(p6,pB,b,val,pp);
276 mgl_dcont_add_pnt(p3,p7,b,val,pp); mgl_dcont_add_pnt(p8,p7,b,val,pp);
277 mgl_dcont_add_pnt(p7,pD,b,val,pp); mgl_dcont_add_pnt(p7,pC,b,val,pp);
278 mgl_dcont_add_pnt(p8,pA,b,val,pp); mgl_dcont_add_pnt(pA,pB,b,val,pp);
279 mgl_dcont_add_pnt(pA,pC,b,val,pp); mgl_dcont_add_pnt(p3,pC,b,val,pp);
280 mgl_dcont_add_pnt(p8,pC,b,val,pp); mgl_dcont_add_pnt(pD,pC,b,val,pp);
281 mgl_dcont_add_pnt(pD,pB,b,val,pp); mgl_dcont_add_pnt(p3,pD,b,val,pp);
282 // now connect points
283 if(pp.size()<2) continue;
284 else if(pp.size()==2)
285 { mglSegment line; line.set(pp[0],pp[1],x,y,z,nboth); lines.push_back(line); }
286 else
287 {
288 mglSegment line;
289 size_t n = pp.size(); // TODO sort by closest in future
290 for(size_t ii=0;ii<n-1;ii++)
291 { line.set(pp[ii],pp[ii+1],x,y,z,nboth); lines.push_back(line); }
292 }
293 }
294 return lines;
295 }
296 //-----------------------------------------------------------------------------
mgl_dcont_gen(HMGL gr,double val,HCDT x,HCDT y,HCDT z,HCDT a,HCDT b,const char * sch,const char * opt)297 void MGL_EXPORT mgl_dcont_gen(HMGL gr, double val, HCDT x, HCDT y, HCDT z, HCDT a, HCDT b, const char *sch, const char *opt)
298 {
299 if(mgl_check_dim3(gr,!mgl_isnboth(x,y,z,a),x,y,z,a,b,"DCont")) return;
300 gr->SaveState(opt);
301 static int cgid=1; gr->StartGroup("DContGen",cgid++);
302
303 int text=0;
304 if(mglchr(sch,'t')) text=1;
305 if(mglchr(sch,'T')) text=2;
306 gr->SetPenPal(sch);
307 mgl_draw_curvs(gr,val,gr->CDef,text,mgl_get_curvs(gr,mgl_get_dlines(val,a,b,x,y,z)));
308 gr->EndGroup();
309 }
310 //-----------------------------------------------------------------------------
mgl_dcont_xyz(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,HCDT a,HCDT b,const char * sch,const char * opt)311 void MGL_EXPORT mgl_dcont_xyz(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, HCDT a, HCDT b, const char *sch, const char *opt)
312 {
313 if(mgl_check_dim3(gr,!mgl_isnboth(x,y,z,a),x,y,z,a,b,"DCont")) return;
314
315 mreal r = gr->SaveState(opt);
316 static int cgid=1; gr->StartGroup("DCont",cgid++);
317
318 int text=0;
319 if(mglchr(sch,'t')) text=1;
320 if(mglchr(sch,'T')) text=2;
321 long s=gr->AddTexture(sch);
322 gr->SetPenPal(sch);
323
324 long Num = mgl_isnan(r)?7:long(r+0.5);
325 if(!v && Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
326 mglData vv(Num);
327 for(long i=0;i<Num;i++) vv.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
328 if(!v) v = &vv;
329
330 #pragma omp parallel for
331 for(long i=0;i<v->GetNx();i++)
332 {
333 mreal val = v->v(i);
334 mgl_draw_curvs(gr,val,gr->GetC(s,val),text, mgl_get_curvs(gr, mgl_get_dlines(val,a,b,x,y,z)));
335 }
336 gr->EndGroup();
337 }
338 //-----------------------------------------------------------------------------
mgl_dcont(HMGL gr,HCDT v,HCDT a,HCDT b,const char * sch,const char * opt)339 void MGL_EXPORT mgl_dcont(HMGL gr, HCDT v, HCDT a, HCDT b, const char *sch, const char *opt)
340 {
341 long n = a->GetNx(), m = a->GetNy(), l = a->GetNz();
342 if(m<2 || n<2 || l<2 || n*m*l!=b->GetNN()) { gr->SetWarn(mglWarnLow,"DCont"); return; }
343 mreal r = gr->SaveState(opt);
344 mglDataV x(n,m,l), y(n,m,l), z(n,m,l);
345 x.Fill(gr->Min.x,gr->Max.x,'x');
346 y.Fill(gr->Min.y,gr->Max.y,'y');
347 z.Fill(gr->Min.z,gr->Max.z,'z');
348 long Num = mgl_isnan(r)?7:long(r+0.5);
349 if(!v && Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
350 mglData vv(Num);
351 for(long i=0;i<Num;i++) vv.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
352 mgl_dcont_xyz(gr,v?v:&vv,&x,&y,&z,a,b,sch,0);
353 }
354 //-----------------------------------------------------------------------------
mgl_dcont_xyz_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,uintptr_t * b,const char * sch,const char * opt,int l,int lo)355 void MGL_EXPORT mgl_dcont_xyz_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, uintptr_t *b, const char *sch, const char *opt,int l,int lo)
356 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
357 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
358 mgl_dcont_xyz(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(z), _DA_(a), _DA_(b), s, o);
359 delete []o; delete []s; }
mgl_dcont_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,uintptr_t * b,const char * sch,const char * opt,int l,int lo)360 void MGL_EXPORT mgl_dcont_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, uintptr_t *b, const char *sch, const char *opt,int l,int lo)
361 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
362 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
363 mgl_dcont(_GR_, _DA_(v), _DA_(a), _DA_(b), s, o);
364 delete []o; delete []s; }
365 //-----------------------------------------------------------------------------
366 //
367 // Cont series
368 //
369 //-----------------------------------------------------------------------------
mgl_get_lines(mreal val,HCDT a,HCDT x,HCDT y,HCDT z,long ak)370 std::vector<mglSegment> MGL_EXPORT mgl_get_lines(mreal val, HCDT a, HCDT x, HCDT y, HCDT z, long ak)
371 {
372 long n=a->GetNx(), m=a->GetNy();
373 std::vector<mglSegment> lines;
374 // first add all possible lines
375 for(long j=0;j<m-1;j++) for(long i=0;i<n-1;i++)
376 {
377 mreal v1=a->v(i,j,ak),v2=a->v(i+1,j,ak),v3=a->v(i,j+1,ak),v4=a->v(i+1,j+1,ak);
378 mreal dl=mgl_d(val,v1,v3),dr=mgl_d(val,v2,v4),dp=mgl_d(val,v1,v2),dn=mgl_d(val,v3,v4);
379 bool added=false;
380 if(v1>val || v4>val)
381 {
382 mglSegment line;
383 if(line.set(0,dl,dn,1,i,j,ak,x,y,z)) { lines.push_back(line); added=true; }
384 if(line.set(1,dr,dp,0,i,j,ak,x,y,z)) { lines.push_back(line); added=true; }
385 }
386 else
387 {
388 mglSegment line;
389 if(line.set(0,dl,dp,0,i,j,ak,x,y,z)) { lines.push_back(line); added=true; }
390 if(line.set(1,dr,dn,1,i,j,ak,x,y,z)) { lines.push_back(line); added=true; }
391 }
392 if(!added) // try to add any other variants
393 {
394 mglSegment line;
395 if(line.set(0,dl,1,dr,i,j,ak,x,y,z)) lines.push_back(line);
396 else if(line.set(dp,0,dn,1,i,j,ak,x,y,z)) lines.push_back(line);
397 else if(line.set(0,dl,dn,1,i,j,ak,x,y,z)) lines.push_back(line);
398 else if(line.set(1,dr,dp,0,i,j,ak,x,y,z)) lines.push_back(line);
399 else if(line.set(0,dl,dp,0,i,j,ak,x,y,z)) lines.push_back(line);
400 else if(line.set(1,dr,dn,1,i,j,ak,x,y,z)) lines.push_back(line);
401 }
402 }
403 return lines;
404 }
405 //-----------------------------------------------------------------------------
mgl_get_curvs(const mglPoint & Min,const mglPoint & Max,std::vector<mglSegment> lines)406 std::vector<mglSegment> MGL_EXPORT mgl_get_curvs(const mglPoint &Min, const mglPoint &Max, std::vector<mglSegment> lines)
407 {
408 long n = lines.size(), m = n;
409 const long nsl=(n>0 && n*n>100)?sqrt(double(n)):10;
410 mreal dxsl = nsl/((Max.x-Min.x)*MGL_FEPSILON), x0 = Min.x;
411 mreal dysl = nsl/((Max.y-Min.y)*MGL_FEPSILON), y0 = Min.y;
412 std::vector<long> *xsl, *ysl;
413 xsl = new std::vector<long>[nsl+1];
414 ysl = new std::vector<long>[nsl+1];
415 for(long i=0;i<n;i++) // group lines by position of its x-coor
416 {
417 long i1 = (lines[i].p1.x-x0)*dxsl, i2 = (lines[i].p2.x-x0)*dxsl;
418 if(i1<0) i1=0;
419 if(i1>nsl) i1=nsl;
420 if(i2<0) i2=0;
421 if(i2>nsl) i2=nsl;
422 if(i1==i2 && i1*(i1-nsl)<=0) xsl[i1].push_back(i);
423 else
424 {
425 if(i1*(i1-nsl)<=0) xsl[i1].push_back(i);
426 if(i2*(i2-nsl)<=0) xsl[i2].push_back(i);
427 }
428 i1 = (lines[i].p1.y-y0)*dysl;
429 i2 = (lines[i].p2.y-y0)*dysl;
430 if(i1<0) i1=0;
431 if(i1>nsl) i1=nsl;
432 if(i2<0) i2=0;
433 if(i2>nsl) i2=nsl;
434 if(i1==i2 && i1*(i1-nsl)<=0) ysl[i1].push_back(i);
435 else
436 {
437 if(i1*(i1-nsl)<=0) ysl[i1].push_back(i);
438 if(i2*(i2-nsl)<=0) ysl[i2].push_back(i);
439 }
440 }
441 size_t xm=0,ym=0;
442 for(long i=0;i<=nsl;i++)
443 {
444 if(xm<xsl[i].size()) xm=xsl[i].size();
445 if(ym<ysl[i].size()) ym=ysl[i].size();
446 }
447 std::vector<mglSegment> curvs;
448 char *used = new char[n]; memset(used,0,n);
449 // create curves from lines
450 while(m>0) // NOTE! This algorithm can be *very* slow!!!
451 {
452 mglSegment curv;
453 bool added = false;
454 for(long i=0;i<n;i++) if(!used[i]) // find any first line segment
455 {
456 curv.before(lines[i].p1);
457 curv.after(lines[i].p2);
458 used[i]=1; m--;
459 added=true; break;
460 }
461 while(added && m>0)
462 {
463 added = false;
464 long i1, i2;
465 if(xm<=ym)
466 { i1 = (curv.p1.x-x0)*dxsl; i2 = (curv.p2.x-x0)*dxsl; }
467 else
468 { i1 = (curv.p1.y-y0)*dysl; i2 = (curv.p2.y-y0)*dysl; }
469 if(i1<0) i1=0;
470 if(i1>nsl) i1=nsl;
471 if(i2<0) i2=0;
472 if(i2>nsl) i2=nsl;
473 const std::vector<long> &isl1=(xm<=ym)?xsl[i1]:ysl[i1];
474 for(size_t i=0;i<isl1.size();i++) // first find continuation of first point
475 {
476 long ii = isl1[i];
477 const mglSegment &l=lines[ii];
478 if(used[ii]) continue;
479 if(l.p1==curv.p1) { curv.before(l.p2); used[ii]=1; m--; added=true; break; }
480 else if(l.p2==curv.p1) { curv.before(l.p1); used[ii]=1; m--; added=true; break; }
481 }
482 const std::vector<long> &isl2=(xm<=ym)?xsl[i2]:ysl[i2];
483 if(m>0) for(size_t i=0;i<isl2.size();i++) // now the same for second point
484 {
485 long ii = isl2[i];
486 const mglSegment &l=lines[ii];
487 if(used[ii]) continue;
488 if(l.p1==curv.p2) { curv.after(l.p2); used[ii]=1; m--; added=true; break; }
489 else if(l.p2==curv.p2) { curv.after(l.p1); used[ii]=1; m--; added=true; break; }
490 }
491 }
492 curvs.push_back(curv);
493 }
494 delete []used; delete []xsl; delete []ysl;
495 return curvs;
496 }
mgl_get_curvs(HMGL gr,const std::vector<mglSegment> & lines)497 std::vector<mglSegment> MGL_EXPORT mgl_get_curvs(HMGL gr, const std::vector<mglSegment> &lines)
498 { return mgl_get_curvs(gr->Min, gr->Max, lines); }
499 //-----------------------------------------------------------------------------
mgl_draw_curvs(HMGL gr,mreal val,mreal c,int text,const std::vector<mglSegment> & curvs)500 void MGL_NO_EXPORT mgl_draw_curvs(HMGL gr, mreal val, mreal c, int text, const std::vector<mglSegment> &curvs)
501 {
502 long pc=0;
503 for(size_t i=0;i<curvs.size();i++) pc += curvs[i].pp.size();
504 gr->Reserve(pc);
505 // fill arguments for other functions
506 long *ff = new long[pc], *nn = new long[pc], m=0;
507 for(size_t i=0;i<curvs.size();i++)
508 {
509 const std::list<mglPoint> &pp=curvs[i].pp;
510 for(std::list<mglPoint>::const_iterator it=pp.begin(); it != pp.end(); ++it)
511 {
512 ff[m] = gr->AddPnt(*it, c);
513 nn[m] = m+1; m++;
514 }
515 nn[m-1]=-1;
516 }
517 if(text && pc>1)
518 {
519 wchar_t wcs[64];
520 mglprintf(wcs,64,L"%4.3g",val);
521 mreal del = 2*gr->TextWidth(wcs,"",-0.5);
522 // find width and height of drawing area
523 mreal ar=gr->GetRatio(), w=gr->GetWidth(), h = gr->GetHeight();
524 ar = (ar>1?1/ar:1)*gr->FontFactor();
525 if(del<ar/2) del = ar/2;
526
527 long m=long(2*w/del)+3, n=long(2*h/del)+3;
528 long *oo=new long[n*m];
529 mreal *rr=new mreal[n*m];
530 for(long i=0;i<n*m;i++) { oo[i]=-1; rr[i]=del*del; }
531 int ii1 = (1664525*pc+1013904223)&0xffff, ii2 = (1664525*ii1+1013904223)&0xffff;
532 mreal x0 = (del*ii1)/0xffff, y0 = (del*ii2)/0xffff;
533 for(long k=0;k<pc;k++) // print label several times if possible
534 {
535 if(nn[k]<0 || ff[k]<0) continue;
536 const mglPoint t = gr->GetPntP(ff[k]);
537 mreal tx = t.x+x0, ty = t.y+y0; // quasi-random shift
538 long i = long(tx/del); tx -= i*del;
539 long j = long(ty/del); ty -= j*del;
540 if(i>=0 && i<m && j>=0 && j<n)
541 {
542 tx = tx*tx+ty*ty; i += m*j;
543 if(rr[i]>tx) { rr[i]=tx; oo[i]=k; }
544 }
545 }
546 for(long i=0;i<n*m;i++) if(oo[i]>=0)
547 mgl_string_curve(gr,oo[i],pc,ff,nn,wcs,text==1?"t:C":"T:C",-0.5);
548 delete []oo; delete []rr;
549 }
550 for(long i=0;i<pc;i++) if(nn[i]>=0) gr->line_plot(ff[i], ff[nn[i]]);
551 delete []nn; delete []ff;
552 }
553 //-----------------------------------------------------------------------------
mgl_data_conts(mreal val,HCDT dat)554 HMDT MGL_EXPORT mgl_data_conts(mreal val, HCDT dat)
555 {
556 mglPoint Min(0,0,0), Max(1,1,1);
557 mglDataV x(dat->GetNx(),dat->GetNy(),dat->GetNz(),0,1,'x');
558 mglDataV y(dat->GetNx(),dat->GetNy(),dat->GetNz(),0,1,'y');
559 mglDataV z(dat->GetNx(),dat->GetNy(),dat->GetNz(),0,1,'z');
560 std::vector<mglSegment> curvs = mgl_get_curvs(Min,Max,mgl_get_lines(val,dat,&x,&y,&z,0));
561 long pc=curvs.size(), m=0;
562 if(pc==0) return NULL;
563 for(size_t i=0;i<curvs.size();i++) pc += curvs[i].pp.size();
564 // fill arguments for other functions
565 HMDT res = new mglData(3,pc);
566 for(size_t i=0;i<curvs.size();i++)
567 {
568 const std::list<mglPoint> &pp=curvs[i].pp;
569 for(std::list<mglPoint>::const_iterator it=pp.begin(); it != pp.end(); ++it)
570 {
571 long i0 = 3*m;
572 res->a[i0] = (*it).x;
573 res->a[1+i0] = (*it).y;
574 res->a[2+i0] = (*it).z;
575 m++;
576 }
577 long i0 = 3*m; m++;
578 res->a[i0] = res->a[1+i0] = res->a[2+i0] = NAN;
579 }
580 return res;
581 }
582 //-----------------------------------------------------------------------------
583 // NOTE! All data MUST have the same size! Only first slice is used!
mgl_cont_genI(HMGL gr,mreal val,HCDT a,HCDT x,HCDT y,HCDT z,mreal c,int text,long ak)584 void MGL_EXPORT mgl_cont_genI(HMGL gr, mreal val, HCDT a, HCDT x, HCDT y, HCDT z, mreal c, int text,long ak)
585 {
586 long n=a->GetNx(), m=a->GetNy();
587 if(n<2 || m<2 || x->GetNx()*x->GetNy()!=n*m || y->GetNx()*y->GetNy()!=n*m || z->GetNx()*z->GetNy()!=n*m)
588 { gr->SetWarn(mglWarnDim,"ContGen"); return; }
589
590 mgl_draw_curvs(gr,val,c,text,mgl_get_curvs(gr,mgl_get_lines(val,a,x,y,z,ak)));
591 }
592 //-----------------------------------------------------------------------------
mgl_cont_gen(HMGL gr,double val,HCDT a,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)593 void MGL_EXPORT mgl_cont_gen(HMGL gr, double val, HCDT a, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
594 {
595 if(mgl_check_dim2(gr,x,y,z,a,"ContGen")) return;
596 gr->SaveState(opt);
597 static int cgid=1; gr->StartGroup("ContGen",cgid++);
598
599 int text=0;
600 if(mglchr(sch,'t')) text=1;
601 if(mglchr(sch,'T')) text=2;
602 gr->SetPenPal(sch);
603 mgl_cont_genI(gr,val,a,x,y,z,gr->CDef,text,0);
604 gr->EndGroup();
605 }
606 //-----------------------------------------------------------------------------
mgl_cont_xy_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)607 void MGL_EXPORT mgl_cont_xy_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
608 {
609 long n=z->GetNx(),m=z->GetNy();
610 if(mgl_check_dim2(gr,x,y,z,0,"Cont")) return;
611
612 gr->SaveState(opt);
613 static int cgid=1; gr->StartGroup("Cont",cgid++);
614
615 int text=0;
616 if(mglchr(sch,'t')) text=1;
617 if(mglchr(sch,'T')) text=2;
618 bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
619 long s=gr->AddTexture(sch);
620 gr->SetPenPal(sch);
621
622 mglData xx, yy;
623 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
624 {
625 xx.Create(n, m); yy.Create(n, m);
626 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
627 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
628 for(long j=0;j<m;j++)
629 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
630 x = &xx; y = &yy;
631 }
632 // x, y -- have the same size z
633 #pragma omp parallel for collapse(2)
634 for(long i=0;i<v->GetNx();i++) for(long j=0;j<z->GetNz();j++)
635 {
636 if(gr->NeedStop()) continue;
637 mreal v0 = v->v(i), z0 = fixed ? gr->Min.z : v0;
638 if(z->GetNz()>1)
639 z0 = gr->Min.z+(gr->Max.z-gr->Min.z)*mreal(j)/(z->GetNz()-1);
640 mglDataV zz(n, m); zz.Fill(z0,z0);
641 mgl_cont_genI(gr,v0,z,x,y,&zz,gr->GetC(s,v0),text,j);
642 }
643 gr->EndGroup();
644 }
645 //-----------------------------------------------------------------------------
mgl_cont_val(HMGL gr,HCDT v,HCDT z,const char * sch,const char * opt)646 void MGL_EXPORT mgl_cont_val(HMGL gr, HCDT v, HCDT z, const char *sch, const char *opt)
647 {
648 long n = z->GetNx(), m = z->GetNy();
649 if(m<2 || n<2) { gr->SetWarn(mglWarnLow,"Cont"); return; }
650 gr->SaveState(opt);
651 mglDataV x(n, m), y(n, m);
652 x.Fill(gr->Min.x,gr->Max.x,'x');
653 y.Fill(gr->Min.y,gr->Max.y,'y');
654 mgl_cont_xy_val(gr,v,&x,&y,z,sch,0);
655 }
656 //-----------------------------------------------------------------------------
657 #define norm(x,y) ((x)*(x)+(y)*(y))
mgl_find_saddle_val(HCDT z)658 std::vector<mreal> MGL_NO_EXPORT mgl_find_saddle_val(HCDT z)
659 {
660 long nx=z->GetNx(), ny=z->GetNy();
661 std::vector<mreal> v;
662 for(long i=1;i<nx-1;i++)
663 {
664 mreal dd=z->vthr(i), x1=z->vthr(i+1), x2=z->vthr(i-1), dp=z->vthr(i+nx);
665 if(dd<=x1 && dd<=x2 && dd>=dp) v.push_back(z->vthr(i));
666 if(dd>=x1 && dd>=x2 && dd<=dp) v.push_back(z->vthr(i));
667 long i0 = i+nx*(ny-1);
668 dd=z->vthr(i0); x1=z->vthr(i0+1); x2=z->vthr(i0-1); dp=z->vthr(i0-nx);
669 if(dd<=x1 && dd<=x2 && dd>=dp) v.push_back(z->vthr(i0));
670 if(dd>=x1 && dd>=x2 && dd<=dp) v.push_back(z->vthr(i0));
671 }
672 for(long j=1;j<ny-1;j++)
673 {
674 long i0 = nx*j;
675 mreal dd=z->vthr(i0), dp=z->vthr(i0+1), y1=z->vthr(i0+nx), y2=z->vthr(i0-nx);
676 if(dd<=dp && dd>=y1 && dd>=y2) v.push_back(z->vthr(i0));
677 if(dd>=dp && dd<=y1 && dd<=y2) v.push_back(z->vthr(i0));
678 i0 = nx*j+nx-1;
679 dd=z->vthr(i0); dp=z->vthr(i0-1); y1=z->vthr(i0+nx); y2=z->vthr(i0-nx);
680 if(dd<=dp && dd>=y1 && dd>=y2) v.push_back(z->vthr(i0));
681 if(dd>=dp && dd<=y1 && dd<=y2) v.push_back(z->vthr(i0));
682 }
683 for(long j=1;j<ny-1;j++) for(long i=1;i<nx-1;i++)
684 {
685 long i0 = i+nx*j; bool ok=false;
686 mreal dd=z->vthr(i0),x1=z->vthr(i0+1),x2=z->vthr(i0-1),y1=z->vthr(i0+nx),y2=z->vthr(i0-nx);
687 if(dd<=x1 && dd<=x2 && dd>=y1 && dd>=y2) ok=true;
688 if(dd>=x1 && dd>=x2 && dd<=y1 && dd<=y2) ok=true;
689 x1=z->vthr(i0+1+nx); x2=z->vthr(i0-1-nx);
690 y1=z->vthr(i0-1+nx); y2=z->vthr(i0+1-nx);
691 if(dd<=x1 && dd<=x2 && dd>=y1 && dd>=y2) ok=true;
692 if(dd>=x1 && dd>=x2 && dd<=y1 && dd<=y2) ok=true;
693 if(ok) v.push_back(z->vthr(i0));
694 }
695 return v;
696 }
mgl_cont_xy(HMGL gr,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)697 void MGL_EXPORT mgl_cont_xy(HMGL gr, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
698 {
699 mreal r = gr->SaveState(opt);
700 if(mglchr(sch,'.'))
701 {
702 mglDataS v; v.dat = mgl_find_saddle_val(z);
703 if(v.dat.size()>0)
704 {
705 std::sort( v.dat.begin(), v.dat.end() );
706 v.dat.erase( std::unique( v.dat.begin(), v.dat.end() ), v.dat.end() );
707 mgl_cont_xy_val(gr,&v,x,y,z,sch,0);
708 }
709 else gr->SetWarn(mglWarnCnt,"Cont");
710 }
711 else
712 {
713 long Num = mgl_isnan(r)?7:long(r+0.5);
714 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
715 mglData v(Num);
716 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
717 mgl_cont_xy_val(gr,&v,x,y,z,sch,0);
718 }
719 }
720 //-----------------------------------------------------------------------------
mgl_cont(HMGL gr,HCDT z,const char * sch,const char * opt)721 void MGL_EXPORT mgl_cont(HMGL gr, HCDT z, const char *sch, const char *opt)
722 {
723 mreal r = gr->SaveState(opt);
724 if(mglchr(sch,'.'))
725 {
726 mglDataS v; v.dat = mgl_find_saddle_val(z);
727 if(v.dat.size()>0)
728 {
729 std::sort( v.dat.begin(), v.dat.end() );
730 v.dat.erase( std::unique( v.dat.begin(), v.dat.end() ), v.dat.end() );
731 mgl_cont_val(gr,&v,z,sch,0);
732 }
733 else gr->SetWarn(mglWarnCnt,"Cont");
734 }
735 else
736 {
737 long Num = mgl_isnan(r)?7:long(r+0.5);
738 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
739 mglData v(Num);
740 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
741 mgl_cont_val(gr,&v,z,sch,0);
742 }
743 }
744 //-----------------------------------------------------------------------------
mgl_cont_xy_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)745 void MGL_EXPORT mgl_cont_xy_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
746 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
747 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
748 mgl_cont_xy_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(a), s, o);
749 delete []o; delete []s; }
750 //-----------------------------------------------------------------------------
mgl_cont_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,const char * opt,int l,int lo)751 void MGL_EXPORT mgl_cont_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
752 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
753 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
754 mgl_cont_val(_GR_, _DA_(v), _DA_(a), s, o); delete []o; delete []s; }
755 //-----------------------------------------------------------------------------
mgl_cont_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)756 void MGL_EXPORT mgl_cont_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
757 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
758 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
759 mgl_cont_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), s, o); delete []o; delete []s; }
760 //-----------------------------------------------------------------------------
mgl_cont_(uintptr_t * gr,uintptr_t * a,const char * sch,const char * opt,int l,int lo)761 void MGL_EXPORT mgl_cont_(uintptr_t *gr, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
762 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
763 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
764 mgl_cont(_GR_, _DA_(a), s, o); delete []o; delete []s; }
765 //-----------------------------------------------------------------------------
766 //
767 // ContF series
768 //
769 //-----------------------------------------------------------------------------
mgl_add_pnt(HMGL gr,mreal d,HCDT x,HCDT y,HCDT z,long i1,long j1,long i2,long j2,mreal c,bool edge)770 long static mgl_add_pnt(HMGL gr, mreal d, HCDT x, HCDT y, HCDT z, long i1, long j1, long i2, long j2, mreal c, bool edge)
771 {
772 long res=-1;
773 if(edge || (d>0 && d<1))
774 {
775 mglPoint p(x->v(i1,j1)*(1-d)+x->v(i2,j2)*d,
776 y->v(i1,j1)*(1-d)+y->v(i2,j2)*d,
777 z->v(i1,j1)*(1-d)+z->v(i2,j2)*d);
778 mglPoint u(x->dvx(i1,j1)*(1-d)+x->dvx(i2,j2)*d,
779 y->dvx(i1,j1)*(1-d)+y->dvx(i2,j2)*d,
780 z->dvx(i1,j1)*(1-d)+z->dvx(i2,j2)*d);
781 mglPoint v(x->dvy(i1,j1)*(1-d)+x->dvy(i2,j2)*d,
782 y->dvy(i1,j1)*(1-d)+y->dvy(i2,j2)*d,
783 z->dvy(i1,j1)*(1-d)+z->dvy(i2,j2)*d);
784 res = gr->AddPnt(p,c,u^v);
785 }
786 return res;
787 }
788 //-----------------------------------------------------------------------------
mgl_add_range(HMGL gr,HCDT a,HCDT x,HCDT y,HCDT z,long i1,long j1,long di,long dj,mreal c,long & u1,long & u2,long ak,mreal v1,mreal v2)789 void static mgl_add_range(HMGL gr, HCDT a, HCDT x, HCDT y, HCDT z, long i1, long j1, long di, long dj, mreal c, long &u1, long &u2, long ak, mreal v1, mreal v2)
790 {
791 long i2=i1+di, j2=j1+dj;
792 mreal f1 = a->v(i1,j1,ak), f2 = a->v(i2,j2,ak), d1, d2;
793 d1 = mgl_d(v1,f1,f2);
794 u1 = mgl_add_pnt(gr,d1,x,y,z,i1,j1,i2,j2,c,false);
795 d2 = mgl_d(v2,f1,f2);
796 u2 = mgl_add_pnt(gr,d2,x,y,z,i1,j1,i2,j2,c,false);
797 if(d1>d2) { j2=u1; u1=u2; u2=j2; }
798 if(u1<0) { u1=u2; u2=-1; }
799 }
800 //-----------------------------------------------------------------------------
mgl_add_edges(HMGL gr,HCDT a,HCDT x,HCDT y,HCDT z,long i1,long j1,long di,long dj,mreal c,long & u1,long & u2,long ak,mreal v1,mreal v2)801 void static mgl_add_edges(HMGL gr, HCDT a, HCDT x, HCDT y, HCDT z, long i1, long j1, long di, long dj, mreal c, long &u1, long &u2, long ak, mreal v1, mreal v2)
802 {
803 long i2=i1+di, j2=j1+dj;
804 u1 = u2 = -1;
805 mreal f1 = a->v(i1,j1,ak), f2 = a->v(i2,j2,ak);
806 if(f1<=v2 && f1>=v1) u1 = mgl_add_pnt(gr,0,x,y,z,i1,j1,i2,j2,c,true);
807 if(f2<=v2 && f2>=v1) u2 = mgl_add_pnt(gr,1,x,y,z,i1,j1,i2,j2,c,true);
808 }
809 //-----------------------------------------------------------------------------
mgl_contf_genI(HMGL gr,mreal v1,mreal v2,HCDT a,HCDT x,HCDT y,HCDT z,mreal c,long ak)810 void MGL_EXPORT mgl_contf_genI(HMGL gr, mreal v1, mreal v2, HCDT a, HCDT x, HCDT y, HCDT z, mreal c, long ak)
811 {
812 long n=a->GetNx(), m=a->GetNy();
813 if(n<2 || m<2 || x->GetNx()*x->GetNy()!=n*m || y->GetNx()*y->GetNy()!=n*m || z->GetNx()*z->GetNy()!=n*m)
814 { gr->SetWarn(mglWarnDim,"ContFGen"); return; }
815
816 gr->Reserve(8*n*m);
817 long *kk = new long[4*n], l1,l2, r1,r2, t1,t2, u1,u2, b1,b2, d1,d2, p[8],num;
818 memset(kk,-1,2*n*sizeof(long));
819 for(long i=0;i<n-1;i++) // add intersection points for first line
820 {
821 mgl_add_range(gr,a,x,y,z, i,0,1,0, c,u1,u2, ak,v1,v2);
822 kk[4*i]=u1; kk[4*i+1]=u2;
823 mgl_add_edges(gr,a,x,y,z, i,0,1,0, c,d1,d2, ak,v1,v2);
824 kk[4*i+2]=d1; kk[4*i+3]=d2;
825 }
826 for(long j=1;j<m;j++) // add intersection points
827 {
828 mgl_add_range(gr,a,x,y,z, 0,j-1,0,1, c,r1,r2, ak,v1,v2);
829 for(long i=0;i<n-1;i++)
830 {
831 l1 = r1; l2 = r2; num=0;
832 t1 = kk[4*i]; t2 = kk[4*i+1];
833 b1 = kk[4*i+2]; b2 = kk[4*i+3];
834 // right edge
835 mgl_add_range(gr,a,x,y,z, i+1,j-1,0,1, c,r1,r2, ak,v1,v2);
836 // top edge
837 mgl_add_range(gr,a,x,y,z, i,j,1,0, c,u1,u2, ak,v1,v2);
838 kk[4*i]=u1; kk[4*i+1]=u2;
839 mgl_add_edges(gr,a,x,y,z, i,j,1,0, c,d1,d2, ak,v1,v2);
840 kk[4*i+2]=d1; kk[4*i+3]=d2;
841 // collect points
842 if(b1>=0) p[num++] = b1;
843 if(t1>=0) p[num++] = t1;
844 if(t2>=0) p[num++] = t2;
845 if(b2>=0) p[num++] = b2;
846 if(r1>=0) p[num++] = r1;
847 if(r2>=0) p[num++] = r2;
848 if(d2>=0) p[num++] = d2;
849 if(u2>=0) p[num++] = u2;
850 if(u1>=0) p[num++] = u1;
851 if(d1>=0) p[num++] = d1;
852 if(l2>=0) p[num++] = l2;
853 if(l1>=0) p[num++] = l1;
854
855 // d1 u1 u2 d2
856 // l2 r2
857 // l1 r1
858 // b1 t1 t2 b2
859
860 // draw it
861 bool b1d2 = a->v(i+1,j,ak)>v2 && a->v(i,j-1,ak)>v2;
862 bool b2d1 = a->v(i,j,ak)>v2 && a->v(i+1,j-1,ak)>v2;
863 // mreal vv = a->linearD(i+0.5,j-0.5,ak,0,0,0);
864 // vv = (vv-v1)*(vv-v2);
865 if(num<3) continue;
866 if(num==4) gr->quad_plot(p[0],p[1],p[3],p[2]);
867 else if(num==3) gr->trig_plot(p[0],p[1],p[2]);
868 else if(num==5)
869 {
870 gr->quad_plot(p[0],p[1],p[3],p[2]);
871 gr->trig_plot(p[0],p[3],p[4]);
872 }
873 else if(num==6)
874 {
875 if(b1>=0 && b2>=0)
876 {
877 gr->quad_plot(b1,b2,l1,r1);
878 gr->quad_plot(l1,r1,u1,u2);
879 }
880 else if(d1>=0 && d2>=0)
881 {
882 gr->quad_plot(d1,d2,l1,r1);
883 gr->quad_plot(l1,r1,t1,t2);
884 }
885 else if(b1>=0 && d2>=0)
886 {
887 if(b2d1)
888 { gr->trig_plot(b1,t1,l1); gr->trig_plot(r1,u1,d2); }
889 else
890 { gr->quad_plot(b1,t1,l1,r1); gr->quad_plot(l1,r1,u1,d2); }
891 }
892 else if(d1>=0 && b2>=0)
893 {
894 if(b1d2)
895 { gr->trig_plot(t1,b2,r1); gr->trig_plot(l1,d1,u1); }
896 else
897 { gr->quad_plot(t1,b2,l1,r1); gr->quad_plot(l1,r1,d1,u1); }
898 }
899 else if(b1>=0 && d1>=0)
900 {
901 gr->quad_plot(b1,d1,t1,u1);
902 gr->quad_plot(t1,u1,r1,r2);
903 }
904 else if(d2>=0 && b2>=0)
905 {
906 gr->quad_plot(d2,b2,u1,t1);
907 gr->quad_plot(t1,u1,l1,l2);
908 }
909 }
910 else if(num==7)
911 {
912 if(b1>=0)
913 {
914 gr->trig_plot(b1,l1,t1); gr->quad_plot(r1,r2,u1,u2);
915 if(!b2d1) gr->quad_plot(l1,t1,u1,r1);
916 }
917 else if(b2>=0)
918 {
919 gr->trig_plot(b2,r1,t1); gr->quad_plot(l1,l2,u2,u1);
920 if(!b1d2) gr->quad_plot(r1,t1,u2,l1);
921 }
922 else if(d2>=0)
923 {
924 gr->trig_plot(d2,r1,u1); gr->quad_plot(l1,l2,t1,t2);
925 if(!b2d1) gr->quad_plot(r1,u1,t2,l2);
926 }
927 else if(d1>=0)
928 {
929 gr->trig_plot(d1,l1,u1); gr->quad_plot(r1,r2,t2,t1);
930 if(!b1d2) gr->quad_plot(l1,u1,t1,r2);
931 }
932 }
933 else if(num==8)
934 {
935 if(b2d1)
936 { if(l2<0) { l2=l1; l1=b1; }
937 if(r2<0) r2=d2;
938 if(t2<0) { t2=t1; t1=b1; }
939 if(u2<0) u2=d2;
940 gr->quad_plot(r1,r2,u1,u2); gr->quad_plot(l1,l2,t1,t2); }
941 else
942 { if(l2<0) l2=d1;
943 if(r2<0) { r2=r1; r1=b2; }
944 if(t2<0) t2=b2;
945 if(u2<0) { u2=u1; u1=d1; }
946 gr->quad_plot(r1,r2,t2,t1); gr->quad_plot(l1,l2,u2,u1); }
947 }
948 }
949 }
950 delete []kk;
951 }
952 //-----------------------------------------------------------------------------
mgl_contf_gen(HMGL gr,double v1,double v2,HCDT a,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)953 void MGL_EXPORT mgl_contf_gen(HMGL gr, double v1, double v2, HCDT a, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
954 {
955 if(mgl_check_dim2(gr,x,y,z,a,"ContFGen")) return;
956 gr->SaveState(opt);
957 static int cgid=1; gr->StartGroup("ContFGen",cgid++);
958
959 gr->SetPenPal(sch);
960 mgl_contf_genI(gr,v1,v2,a,x,y,z,gr->CDef,0);
961 gr->EndGroup();
962 }
963 //-----------------------------------------------------------------------------
mgl_contf_xy_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)964 void MGL_EXPORT mgl_contf_xy_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
965 {
966 long n=z->GetNx(),m=z->GetNy();
967 if(mgl_check_dim2(gr,x,y,z,0,"ContF")) return;
968
969 gr->SaveState(opt);
970 static int cgid=1; gr->StartGroup("ContF",cgid++);
971 long s=gr->AddTexture(sch);
972
973 bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
974 mglData xx, yy;
975 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
976 {
977 xx.Create(n, m); yy.Create(n, m);
978 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
979 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
980 for(long j=0;j<m;j++)
981 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
982 x = &xx; y = &yy;
983 }
984 // x, y -- have the same size z
985 #pragma omp parallel for collapse(2)
986 for(long i=0;i<v->GetNx()-1;i++) for(long j=0;j<z->GetNz();j++)
987 {
988 if(gr->NeedStop()) continue;
989 mreal v0 = v->v(i), z0 = fixed ? gr->Min.z : v0;
990 if(z->GetNz()>1)
991 z0 = gr->Min.z+(gr->Max.z-gr->Min.z)*mreal(j)/(z->GetNz()-1);
992 mglDataV zz(n, m); zz.Fill(z0,z0);
993 mgl_contf_genI(gr,v0,v->v(i+1),z,x,y,&zz,gr->GetC(s,v0),j);
994 }
995 gr->EndGroup();
996 }
997 //-----------------------------------------------------------------------------
mgl_contf_val(HMGL gr,HCDT v,HCDT z,const char * sch,const char * opt)998 void MGL_EXPORT mgl_contf_val(HMGL gr, HCDT v, HCDT z, const char *sch, const char *opt)
999 {
1000 long n = z->GetNx(), m = z->GetNy();
1001 if(n<2 || m<2) { gr->SetWarn(mglWarnLow,"Cont"); return; }
1002 gr->SaveState(opt);
1003 mglDataV x(n, m), y(n, m);
1004 x.Fill(gr->Min.x,gr->Max.x,'x');
1005 y.Fill(gr->Min.y,gr->Max.y,'y');
1006 mgl_contf_xy_val(gr,v,&x,&y,z,sch,0);
1007 }
1008 //-----------------------------------------------------------------------------
mgl_contf_xy(HMGL gr,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1009 void MGL_EXPORT mgl_contf_xy(HMGL gr, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1010 {
1011 mreal r = gr->SaveState(opt);
1012 long Num = mgl_isnan(r)?7:long(r+0.5);
1013 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
1014 mglDataV v(Num+2); v.Fill(gr->Min.c, gr->Max.c);
1015 mgl_contf_xy_val(gr,&v,x,y,z,sch,0);
1016 }
1017 //-----------------------------------------------------------------------------
mgl_contf(HMGL gr,HCDT z,const char * sch,const char * opt)1018 void MGL_EXPORT mgl_contf(HMGL gr, HCDT z, const char *sch, const char *opt)
1019 {
1020 mreal r = gr->SaveState(opt);
1021 long Num = mgl_isnan(r)?7:long(r+0.5);
1022 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
1023 mglDataV v(Num+2); v.Fill(gr->Min.c, gr->Max.c);
1024 mgl_contf_val(gr,&v,z,sch,0);
1025 }
1026 //-----------------------------------------------------------------------------
mgl_contf_xy_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1027 void MGL_EXPORT mgl_contf_xy_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1028 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1029 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1030 mgl_contf_xy_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(a), s, o);
1031 delete []o; delete []s; }
1032 //-----------------------------------------------------------------------------
mgl_contf_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1033 void MGL_EXPORT mgl_contf_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1034 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1035 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1036 mgl_contf_val(_GR_, _DA_(v), _DA_(a), s, o); delete []o; delete []s; }
1037 //-----------------------------------------------------------------------------
mgl_contf_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1038 void MGL_EXPORT mgl_contf_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1039 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1040 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1041 mgl_contf_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), s, o);
1042 delete []o; delete []s; }
1043 //-----------------------------------------------------------------------------
mgl_contf_(uintptr_t * gr,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1044 void MGL_EXPORT mgl_contf_(uintptr_t *gr, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1045 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1046 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1047 mgl_contf(_GR_, _DA_(a), s, o); delete []o; delete []s; }
1048 //-----------------------------------------------------------------------------
1049 //
1050 // ContP series
1051 //
1052 //-----------------------------------------------------------------------------
mgl_contp_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,const char * opt)1053 void MGL_EXPORT mgl_contp_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, const char *opt)
1054 {
1055 long n=z->GetNx(),m=z->GetNy();
1056 if(mgl_check_dim2(gr,x,y,z,a,"Cont")) return;
1057
1058 gr->SaveState(opt);
1059 static int cgid=1; gr->StartGroup("Cont",cgid++);
1060
1061 int text=0;
1062 if(mglchr(sch,'t')) text=1;
1063 if(mglchr(sch,'T')) text=2;
1064 bool fill = mglchr(sch,'f');
1065 long s=gr->AddTexture(sch);
1066 gr->SetPenPal(sch);
1067
1068 mglData xx, yy;
1069 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
1070 {
1071 xx.Create(n, m); yy.Create(n, m);
1072 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
1073 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
1074 for(long j=0;j<m;j++)
1075 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
1076 x = &xx; y = &yy;
1077 }
1078 // x, y -- have the same size z
1079 #pragma omp parallel for collapse(2)
1080 for(long i=0;i<v->GetNx();i++) for(long j=0;j<a->GetNz();j++)
1081 {
1082 if(gr->NeedStop()) continue;
1083 if(fill)
1084 mgl_contf_genI(gr,v->v(i),v->v(i+1),a,x,y,z,gr->GetC(s,v->v(i)),j);
1085 else
1086 mgl_cont_genI(gr,v->v(i),a,x,y,z,gr->GetC(s,v->v(i)),text,j);
1087 }
1088 gr->EndGroup();
1089 }
1090 //-----------------------------------------------------------------------------
mgl_contp(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,const char * opt)1091 void MGL_EXPORT mgl_contp(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, const char *opt)
1092 {
1093 mreal r = gr->SaveState(opt);
1094 long Num = mgl_isnan(r)?7:long(r+0.5);
1095 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
1096 mglData v(Num);
1097 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1098 mgl_contp_val(gr,&v,x,y,z,a,sch,0);
1099 }
1100 //-----------------------------------------------------------------------------
mgl_contp_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1101 void MGL_EXPORT mgl_contp_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1102 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1103 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1104 mgl_contp_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, o);
1105 delete []o; delete []s; }
1106 //-----------------------------------------------------------------------------
mgl_contp_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1107 void MGL_EXPORT mgl_contp_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1108 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1109 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1110 mgl_contp(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, o);
1111 delete []o; delete []s; }
1112 //-----------------------------------------------------------------------------
1113 //
1114 // ContD series
1115 //
1116 //-----------------------------------------------------------------------------
mgl_get_ncol(const char * sch,char * res)1117 int static mgl_get_ncol(const char *sch, char *res)
1118 {
1119 long j=0;
1120 if(sch) for(long i=0;sch[i]&&sch[i]!=':';i++) if(strchr(MGL_COLORS,sch[i]))
1121 { if(res) res[j]=sch[i]; j++; }
1122 return j?j:strlen(MGL_DEF_PAL);
1123 }
1124 //-----------------------------------------------------------------------------
mgl_contd_xy_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1125 void MGL_EXPORT mgl_contd_xy_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1126 {
1127 long j=0,n=z->GetNx(),m=z->GetNy();
1128 if(mgl_check_dim2(gr,x,y,z,0,"ContD")) return;
1129
1130 gr->SaveState(opt);
1131 static int cgid=1; gr->StartGroup("ContD",cgid++);
1132
1133 bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
1134 if(sch) for(long i=0;sch[i];i++) if(strchr(MGL_COLORS,sch[i])) j++;
1135 if(j==0) sch = MGL_DEF_PAL;
1136 long s = gr->AddTexture(sch,1);
1137 int nc = gr->GetNumPal(s*256);
1138 mglData xx, yy;
1139 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
1140 {
1141 xx.Create(n, m); yy.Create(n, m);
1142 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
1143 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
1144 for(long j=0;j<m;j++)
1145 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
1146 x = &xx; y = &yy;
1147 }
1148 // x, y -- have the same size z
1149 mreal dc = nc>1 ? 1/(MGL_FEPSILON*(nc-1)) : 0;
1150 #pragma omp parallel for collapse(2)
1151 for(long i=0;i<v->GetNx()-1;i++) for(long j=0;j<z->GetNz();j++)
1152 {
1153 if(gr->NeedStop()) continue;
1154 mreal v0 = v->v(i), z0 = fixed ? gr->Min.z : v0;
1155 if(z->GetNz()>1)
1156 z0 = gr->Min.z+(gr->Max.z-gr->Min.z)*mreal(j)/(z->GetNz()-1);
1157 mglDataV zz(n, m); zz.Fill(z0,z0);
1158 mgl_contf_genI(gr,v0,v->v(i+1),z,x,y,&zz,s+(i%nc)*dc,j);
1159 }
1160 gr->EndGroup();
1161 }
1162 //-----------------------------------------------------------------------------
mgl_contd_val(HMGL gr,HCDT v,HCDT z,const char * sch,const char * opt)1163 void MGL_EXPORT mgl_contd_val(HMGL gr, HCDT v, HCDT z, const char *sch, const char *opt)
1164 {
1165 long n = z->GetNx(), m = z->GetNy();
1166 if(n<2 || m<2) { gr->SetWarn(mglWarnLow,"ContD"); return; }
1167 gr->SaveState(opt);
1168 mglDataV x(n, m), y(n, m);
1169 x.Fill(gr->Min.x,gr->Max.x,'x');
1170 y.Fill(gr->Min.y,gr->Max.y,'y');
1171 mgl_contd_xy_val(gr,v,&x,&y,z,sch,0);
1172 }
1173 //-----------------------------------------------------------------------------
mgl_contd_xy(HMGL gr,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1174 void MGL_EXPORT mgl_contd_xy(HMGL gr, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1175 {
1176 gr->SaveState(opt);
1177 mglDataV v(mgl_get_ncol(sch,0)+1);
1178 v.Fill(gr->Min.c, gr->Max.c);
1179 mgl_contd_xy_val(gr,&v,x,y,z,sch,0);
1180 }
1181 //-----------------------------------------------------------------------------
mgl_contd(HMGL gr,HCDT z,const char * sch,const char * opt)1182 void MGL_EXPORT mgl_contd(HMGL gr, HCDT z, const char *sch, const char *opt)
1183 {
1184 gr->SaveState(opt);
1185 mglDataV v(mgl_get_ncol(sch,0)+1);
1186 v.Fill(gr->Min.c, gr->Max.c);
1187 mgl_contd_val(gr,&v,z,sch,0);
1188 }
1189 //-----------------------------------------------------------------------------
mgl_contd_xy_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1190 void MGL_EXPORT mgl_contd_xy_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1191 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1192 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1193 mgl_contd_xy_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(a), s, o);
1194 delete []o; delete []s; }
1195 //-----------------------------------------------------------------------------
mgl_contd_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1196 void MGL_EXPORT mgl_contd_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1197 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1198 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1199 mgl_contd_val(_GR_, _DA_(v), _DA_(a), s, o); delete []o; delete []s; }
1200 //-----------------------------------------------------------------------------
mgl_contd_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1201 void MGL_EXPORT mgl_contd_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1202 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1203 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1204 mgl_contd_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), s, o); delete []o; delete []s; }
1205 //-----------------------------------------------------------------------------
mgl_contd_(uintptr_t * gr,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1206 void MGL_EXPORT mgl_contd_(uintptr_t *gr, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1207 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1208 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1209 mgl_contd(_GR_, _DA_(a), s, o); delete []o; delete []s; }
1210 //-----------------------------------------------------------------------------
1211 //
1212 // ContV series
1213 //
1214 //-----------------------------------------------------------------------------
1215 // NOTE! All data MUST have the same size! Only first slice is used!
mgl_contv_genI(HMGL gr,mreal val,mreal dval,HCDT a,HCDT x,HCDT y,HCDT z,mreal c,long ak)1216 void MGL_EXPORT mgl_contv_genI(HMGL gr, mreal val, mreal dval, HCDT a, HCDT x, HCDT y, HCDT z, mreal c, long ak)
1217 {
1218 long n=a->GetNx(), m=a->GetNy();
1219 if(n<2 || m<2 || x->GetNx()*x->GetNy()!=n*m || y->GetNx()*y->GetNy()!=n*m || z->GetNx()*z->GetNy()!=n*m)
1220 { gr->SetWarn(mglWarnDim,"ContGen"); return; }
1221
1222 const std::vector<mglSegment> curvs = mgl_get_curvs(gr,mgl_get_lines(val,a,x,y,z,ak));
1223 for(size_t i=0;i<curvs.size();i++)
1224 {
1225 const std::list<mglPoint> &pp=curvs[i].pp;
1226 long f2=-1,g2=-1;
1227 for(std::list<mglPoint>::const_iterator it=pp.begin(); it != pp.end(); ++it)
1228 {
1229 mglPoint p=*it,q(p.y,-p.x);
1230 long f1 = f2; f2 = gr->AddPnt(p,c,q); p.z+=dval;
1231 long g1 = g2; g2 = gr->AddPnt(p,c,q);
1232 gr->quad_plot(f1,g1,f2,g2);
1233 }
1234 }
1235 }
1236 //-----------------------------------------------------------------------------
mgl_contv_xy_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1237 void MGL_EXPORT mgl_contv_xy_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1238 {
1239 long n=z->GetNx(),m=z->GetNy();
1240 if(mgl_check_dim2(gr,x,y,z,0,"ContV")) return;
1241
1242 gr->SaveState(opt);
1243 static int cgid=1; gr->StartGroup("ContV",cgid++);
1244 bool fixed=(mglchr(sch,'_')) || (gr->Min.z==gr->Max.z);
1245 long s=gr->AddTexture(sch);
1246 gr->SetPenPal(sch);
1247
1248 mglData xx, yy;
1249 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
1250 {
1251 xx.Create(n, m); yy.Create(n, m);
1252 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
1253 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
1254 for(long j=0;j<m;j++)
1255 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
1256 x = &xx; y = &yy;
1257 }
1258 // x, y -- have the same size z
1259 #pragma omp parallel for collapse(2)
1260 for(long i=0;i<v->GetNx();i++) for(long j=0;j<z->GetNz();j++)
1261 {
1262 if(gr->NeedStop()) continue;
1263 mreal v0 = v->v(i), z0 = fixed ? gr->Min.z : v0;
1264 if(z->GetNz()>1) z0 = gr->Min.z+(gr->Max.z-gr->Min.z)*mreal(j)/(z->GetNz()-1);
1265 mglDataV zz(n, m); zz.Fill(z0,z0);
1266 mreal dv = (gr->Max.c-gr->Min.c)/8;
1267 if(i>0) dv = v->v(i-1)-v->v(i);
1268 else if(i<v->GetNx()-1) dv = v->v(i)-v->v(i+1);
1269 if(fixed) dv=-dv;
1270 mgl_contv_genI(gr,v0,dv,z,x,y,&zz,gr->GetC(s,v0),j);
1271 }
1272 gr->EndGroup();
1273 }
1274 //-----------------------------------------------------------------------------
mgl_contv_val(HMGL gr,HCDT v,HCDT z,const char * sch,const char * opt)1275 void MGL_EXPORT mgl_contv_val(HMGL gr, HCDT v, HCDT z, const char *sch, const char *opt)
1276 {
1277 long n = z->GetNx(), m = z->GetNy();
1278 if(n<2 || m<2) { gr->SetWarn(mglWarnLow,"Cont"); return; }
1279 gr->SaveState(opt);
1280 mglDataV x(n, m), y(n, m);
1281 x.Fill(gr->Min.x,gr->Max.x,'x');
1282 y.Fill(gr->Min.y,gr->Max.y,'y');
1283 mgl_contv_xy_val(gr,v,&x,&y,z,sch,0);
1284 }
1285 //-----------------------------------------------------------------------------
mgl_contv_xy(HMGL gr,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1286 void MGL_EXPORT mgl_contv_xy(HMGL gr, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1287 {
1288 mreal r = gr->SaveState(opt);
1289 long Num = mgl_isnan(r)?7:long(r+0.5);
1290 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
1291 mglData v(Num);
1292 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1293 mgl_contv_xy_val(gr,&v,x,y,z,sch,0);
1294 }
1295 //-----------------------------------------------------------------------------
mgl_contv(HMGL gr,HCDT z,const char * sch,const char * opt)1296 void MGL_EXPORT mgl_contv(HMGL gr, HCDT z, const char *sch, const char *opt)
1297 {
1298 mreal r = gr->SaveState(opt);
1299 long Num = mgl_isnan(r)?7:long(r+0.5);
1300 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont"); return; }
1301 mglData v(Num);
1302 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1303 mgl_contv_val(gr,&v,z,sch,0);
1304 }
1305 //-----------------------------------------------------------------------------
mgl_contv_xy_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1306 void MGL_EXPORT mgl_contv_xy_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1307 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1308 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1309 mgl_contv_xy_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(a), s, o);
1310 delete []o; delete []s; }
1311 //-----------------------------------------------------------------------------
mgl_contv_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1312 void MGL_EXPORT mgl_contv_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1313 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1314 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1315 mgl_contv_val(_GR_, _DA_(v), _DA_(a), s, o); delete []o; delete []s; }
1316 //-----------------------------------------------------------------------------
mgl_contv_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1317 void MGL_EXPORT mgl_contv_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1318 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1319 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1320 mgl_contv_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), s, o); delete []o; delete []s; }
1321 //-----------------------------------------------------------------------------
mgl_contv_(uintptr_t * gr,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1322 void MGL_EXPORT mgl_contv_(uintptr_t *gr, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1323 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1324 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1325 mgl_contv(_GR_, _DA_(a), s, o); delete []o; delete []s; }
1326 //-----------------------------------------------------------------------------
1327 //
1328 // Cont3 series
1329 //
1330 //-----------------------------------------------------------------------------
1331 struct _mgl_slice { mglData x,y,z,a; };
1332 //-----------------------------------------------------------------------------
mgl_get_slice(_mgl_slice & s,HCDT x,HCDT y,HCDT z,HCDT a,char dir,mreal d,bool both)1333 void MGL_NO_EXPORT mgl_get_slice(_mgl_slice &s, HCDT x, HCDT y, HCDT z, HCDT a, char dir, mreal d, bool both)
1334 {
1335 long n=a->GetNx(),m=a->GetNy(),l=a->GetNz(), nx=1,ny=1,p;
1336
1337 if(dir=='x') { nx = m; ny = l; if(d<0) d = n/2.; }
1338 if(dir=='y') { nx = n; ny = l; if(d<0) d = m/2.; }
1339 if(dir=='z') { nx = n; ny = m; if(d<0) d = l/2.; }
1340 s.x.Create(nx,ny); s.y.Create(nx,ny);
1341 s.z.Create(nx,ny); s.a.Create(nx,ny);
1342 p = long(d); d -= p;
1343 if(dir=='x' && p>=n-1) { d+=p-n+2; p=n-2; }
1344 if(dir=='y' && p>=m-1) { d+=p-m+2.; p=m-2; }
1345 if(dir=='z' && p>=l-1) { d+=p-l+2; p=l-2; }
1346 mreal v;
1347
1348 if(both)
1349 {
1350 if(dir=='x')
1351 #pragma omp parallel for collapse(2)
1352 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1353 {
1354 long i0 = i+nx*j;
1355 s.x.a[i0] = x->v(p,i,j)*(1-d) + x->v(p+1,i,j)*d;
1356 s.y.a[i0] = y->v(p,i,j)*(1-d) + y->v(p+1,i,j)*d;
1357 s.z.a[i0] = z->v(p,i,j)*(1-d) + z->v(p+1,i,j)*d;
1358 s.a.a[i0] = a->v(p,i,j)*(1-d) + a->v(p+1,i,j)*d;
1359 }
1360 if(dir=='y')
1361 #pragma omp parallel for collapse(2)
1362 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1363 {
1364 long i0 = i+nx*j;
1365 s.x.a[i0] = x->v(i,p,j)*(1-d) + x->v(i,p+1,j)*d;
1366 s.y.a[i0] = y->v(i,p,j)*(1-d) + y->v(i,p+1,j)*d;
1367 s.z.a[i0] = z->v(i,p,j)*(1-d) + z->v(i,p+1,j)*d;
1368 s.a.a[i0] = a->v(i,p,j)*(1-d) + a->v(i,p+1,j)*d;
1369 }
1370 if(dir=='z')
1371 #pragma omp parallel for collapse(2)
1372 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1373 {
1374 long i0 = i+nx*j;
1375 s.x.a[i0] = x->v(i,j,p)*(1-d) + x->v(i,j,p+1)*d;
1376 s.y.a[i0] = y->v(i,j,p)*(1-d) + y->v(i,j,p+1)*d;
1377 s.z.a[i0] = z->v(i,j,p)*(1-d) + z->v(i,j,p+1)*d;
1378 s.a.a[i0] = a->v(i,j,p)*(1-d) + a->v(i,j,p+1)*d;
1379 }
1380 }
1381 else // x, y, z -- vectors
1382 {
1383 if(dir=='x')
1384 {
1385 v = x->v(p)*(1-d)+x->v(p+1)*d;
1386 #pragma omp parallel for collapse(2)
1387 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1388 {
1389 long i0 = i+nx*j; s.x.a[i0] = v;
1390 s.y.a[i0] = y->v(i); s.z.a[i0] = z->v(j);
1391 s.a.a[i0] = a->v(p,i,j)*(1-d) + a->v(p+1,i,j)*d;
1392 }
1393 }
1394 if(dir=='y')
1395 {
1396 v = y->v(p)*(1-d)+y->v(p+1)*d;
1397 #pragma omp parallel for collapse(2)
1398 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1399 {
1400 long i0 = i+nx*j; s.y.a[i0] = v;
1401 s.x.a[i0] = x->v(i); s.z.a[i0] = z->v(j);
1402 s.a.a[i0] = a->v(i,p,j)*(1-d) + a->v(i,p+1,j)*d;
1403 }
1404 }
1405 if(dir=='z')
1406 {
1407 v = z->v(p)*(1-d)+z->v(p+1)*d;
1408 #pragma omp parallel for collapse(2)
1409 for(long j=0;j<ny;j++) for(long i=0;i<nx;i++)
1410 {
1411 long i0 = i+nx*j; s.z.a[i0] = v;
1412 s.x.a[i0] = x->v(i); s.y.a[i0] = y->v(j);
1413 s.a.a[i0] = a->v(i,j,p)*(1-d) + a->v(i,j,p+1)*d;
1414 }
1415 }
1416 }
1417 }
1418 //-----------------------------------------------------------------------------
mgl_cont3_xyz_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1419 void MGL_EXPORT mgl_cont3_xyz_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1420 {
1421 bool both = mgl_isboth(x,y,z,a);
1422 if(mgl_check_dim3(gr,both,x,y,z,a,0,"Cont3")) return;
1423
1424 gr->SaveState(opt);
1425 static int cgid=1; gr->StartGroup("Cont3",cgid++);
1426 char dir='y';
1427 if(mglchr(sch,'x')) dir='x';
1428 if(mglchr(sch,'z')) dir='z';
1429
1430 int text=0;
1431 if(mglchr(sch,'t')) text=1;
1432 if(mglchr(sch,'T')) text=2;
1433 long ss=gr->AddTexture(sch);
1434 gr->SetPenPal(sch);
1435
1436 _mgl_slice s;
1437 mgl_get_slice(s,x,y,z,a,dir,sVal,both);
1438 #pragma omp parallel for
1439 for(long i=0;i<v->GetNx();i++)
1440 {
1441 mreal v0 = v->v(i);
1442 mgl_cont_genI(gr,v0,&s.a,&s.x,&s.y,&s.z,gr->GetC(ss,v0),text,0);
1443 }
1444 gr->EndGroup();
1445 }
1446 //-----------------------------------------------------------------------------
mgl_cont3_val(HMGL gr,HCDT v,HCDT a,const char * sch,double sVal,const char * opt)1447 void MGL_EXPORT mgl_cont3_val(HMGL gr, HCDT v, HCDT a, const char *sch, double sVal, const char *opt)
1448 {
1449 gr->SaveState(opt);
1450 mglDataV x(a->GetNx()), y(a->GetNy()),z(a->GetNz());
1451 x.Fill(gr->Min.x,gr->Max.x);
1452 y.Fill(gr->Min.y,gr->Max.y);
1453 z.Fill(gr->Min.z,gr->Max.z);
1454 mgl_cont3_xyz_val(gr,v,&x,&y,&z,a,sch,sVal,0);
1455 }
1456 //-----------------------------------------------------------------------------
mgl_cont3_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1457 void MGL_EXPORT mgl_cont3_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1458 {
1459 mreal r = gr->SaveState(opt);
1460 long Num = mgl_isnan(r)?7:long(r+0.5);
1461 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont3"); return; }
1462 mglData v(Num);
1463 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1464 mgl_cont3_xyz_val(gr,&v,x,y,z,a,sch,sVal,0);
1465 }
1466 //-----------------------------------------------------------------------------
mgl_cont3(HMGL gr,HCDT a,const char * sch,double sVal,const char * opt)1467 void MGL_EXPORT mgl_cont3(HMGL gr, HCDT a, const char *sch, double sVal, const char *opt)
1468 {
1469 mreal r = gr->SaveState(opt);
1470 long Num = mgl_isnan(r)?7:long(r+0.5);
1471 if(Num<1) { gr->SetWarn(mglWarnCnt,"Cont3"); return; }
1472 mglData v(Num);
1473 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1474 mgl_cont3_val(gr,&v,a,sch,sVal,0);
1475 }
1476 //-----------------------------------------------------------------------------
mgl_cont3_xyz_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1477 void MGL_EXPORT mgl_cont3_xyz_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1478 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1479 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1480 mgl_cont3_xyz_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1481 delete []o; delete []s; }
1482 //-----------------------------------------------------------------------------
mgl_cont3_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1483 void MGL_EXPORT mgl_cont3_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1484 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1485 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1486 mgl_cont3_val(_GR_, _DA_(v), _DA_(a), s, *sVal, o); delete []o; delete []s; }
1487 //-----------------------------------------------------------------------------
mgl_cont3_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1488 void MGL_EXPORT mgl_cont3_xyz_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1489 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1490 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1491 mgl_cont3_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1492 delete []o; delete []s; }
1493 //-----------------------------------------------------------------------------
mgl_cont3_(uintptr_t * gr,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1494 void MGL_EXPORT mgl_cont3_(uintptr_t *gr, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1495 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1496 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1497 mgl_cont3(_GR_, _DA_(a), s, *sVal, o); delete []o; delete []s; }
1498 //-----------------------------------------------------------------------------
1499 //
1500 // Dens3 series
1501 //
1502 //-----------------------------------------------------------------------------
1503 void MGL_NO_EXPORT mgl_surf_gen(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT c, HCDT a, const char *sch);
mgl_dens3_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1504 void MGL_EXPORT mgl_dens3_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1505 {
1506 bool both = mgl_isboth(x,y,z,a);
1507 if(mgl_check_dim3(gr,both,x,y,z,a,0,"Dens3")) return;
1508
1509 gr->SaveState(opt);
1510 static int cgid=1; gr->StartGroup("Dens3",cgid++);
1511 char dir='y';
1512 if(mglchr(sch,'x')) dir='x';
1513 if(mglchr(sch,'z')) dir='z';
1514
1515 _mgl_slice s;
1516 mgl_get_slice(s,x,y,z,a,dir,sVal,both);
1517 mgl_surf_gen(gr, &s.x,&s.y,&s.z,&s.a, 0, sch);
1518 // mgl_surfc_xy(gr,&s.x,&s.y,&s.z,&s.a,sch,0);
1519 gr->EndGroup();
1520 }
1521 //-----------------------------------------------------------------------------
mgl_dens3(HMGL gr,HCDT a,const char * sch,double sVal,const char * opt)1522 void MGL_EXPORT mgl_dens3(HMGL gr, HCDT a, const char *sch, double sVal, const char *opt)
1523 {
1524 gr->SaveState(opt);
1525 mglDataV x(a->GetNx()), y(a->GetNy()),z(a->GetNz());
1526 x.Fill(gr->Min.x,gr->Max.x);
1527 y.Fill(gr->Min.y,gr->Max.y);
1528 z.Fill(gr->Min.z,gr->Max.z);
1529 mgl_dens3_xyz(gr,&x,&y,&z,a,sch,sVal,0);
1530 }
1531 //-----------------------------------------------------------------------------
mgl_dens3_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1532 void MGL_EXPORT mgl_dens3_xyz_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1533 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1534 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1535 mgl_dens3_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1536 delete []o; delete []s; }
1537 //-----------------------------------------------------------------------------
mgl_dens3_(uintptr_t * gr,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1538 void MGL_EXPORT mgl_dens3_(uintptr_t *gr, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1539 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1540 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1541 mgl_dens3(_GR_, _DA_(a), s, *sVal, o); delete []o; delete []s; }
1542 //-----------------------------------------------------------------------------
1543 //
1544 // Grid3 series
1545 //
1546 //-----------------------------------------------------------------------------
mgl_grid3_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1547 void MGL_EXPORT mgl_grid3_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1548 {
1549 bool both = mgl_isboth(x,y,z,a);
1550 if(mgl_check_dim3(gr,both,x,y,z,a,0,"Grid3")) return;
1551
1552 gr->SaveState(opt);
1553 static int cgid=1; gr->StartGroup("Grid3",cgid++);
1554 char dir='y';
1555 if(mglchr(sch,'x')) dir='x';
1556 if(mglchr(sch,'z')) dir='z';
1557
1558 _mgl_slice s;
1559 mgl_get_slice(s,x,y,z,a,dir,sVal,both);
1560 mgl_mesh_xy(gr,&s.x,&s.y,&s.z,sch,0);
1561 gr->EndGroup();
1562 }
1563 //-----------------------------------------------------------------------------
mgl_grid3(HMGL gr,HCDT a,const char * sch,double sVal,const char * opt)1564 void MGL_EXPORT mgl_grid3(HMGL gr, HCDT a, const char *sch, double sVal, const char *opt)
1565 {
1566 gr->SaveState(opt);
1567 mglDataV x(a->GetNx()), y(a->GetNy()), z(a->GetNz());
1568 x.Fill(gr->Min.x,gr->Max.x);
1569 y.Fill(gr->Min.y,gr->Max.y);
1570 z.Fill(gr->Min.z,gr->Max.z);
1571 mgl_grid3_xyz(gr,&x,&y,&z,a,sch,sVal,0);
1572 }
1573 //-----------------------------------------------------------------------------
mgl_grid3_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1574 void MGL_EXPORT mgl_grid3_xyz_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1575 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1576 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1577 mgl_grid3_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1578 delete []o; delete []s; }
1579 //-----------------------------------------------------------------------------
mgl_grid3_(uintptr_t * gr,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1580 void MGL_EXPORT mgl_grid3_(uintptr_t *gr, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1581 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1582 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1583 mgl_grid3(_GR_, _DA_(a), s, *sVal, o); delete []o; delete []s; }
1584 //-----------------------------------------------------------------------------
1585 //
1586 // ContF3 series
1587 //
1588 //-----------------------------------------------------------------------------
mgl_contf3_xyz_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1589 void MGL_EXPORT mgl_contf3_xyz_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1590 {
1591 bool both = mgl_isboth(x,y,z,a);
1592 if(mgl_check_dim3(gr,both,x,y,z,a,0,"ContF3")) return;
1593
1594 gr->SaveState(opt);
1595 static int cgid=1; gr->StartGroup("ContF3",cgid++);
1596 char dir='y';
1597 if(mglchr(sch,'x')) dir='x';
1598 if(mglchr(sch,'z')) dir='z';
1599
1600 long ss=gr->AddTexture(sch);
1601 _mgl_slice s;
1602 mgl_get_slice(s,x,y,z,a,dir,sVal,both);
1603 #pragma omp parallel for
1604 for(long i=0;i<v->GetNx()-1;i++)
1605 {
1606 mreal v0 = v->v(i);
1607 mgl_contf_genI(gr,v0,v->v(i+1),&s.a,&s.x,&s.y,&s.z,gr->GetC(ss,v0),0);
1608 }
1609 gr->EndGroup();
1610 }
1611 //-----------------------------------------------------------------------------
mgl_contf3_val(HMGL gr,HCDT v,HCDT a,const char * sch,double sVal,const char * opt)1612 void MGL_EXPORT mgl_contf3_val(HMGL gr, HCDT v, HCDT a, const char *sch, double sVal, const char *opt)
1613 {
1614 gr->SaveState(opt);
1615 mglDataV x(a->GetNx()), y(a->GetNy()),z(a->GetNz());
1616 x.Fill(gr->Min.x,gr->Max.x);
1617 y.Fill(gr->Min.y,gr->Max.y);
1618 z.Fill(gr->Min.z,gr->Max.z);
1619 mgl_contf3_xyz_val(gr,v,&x,&y,&z,a,sch,sVal,0);
1620 }
1621 //-----------------------------------------------------------------------------
mgl_contf3_xyz(HMGL gr,HCDT x,HCDT y,HCDT z,HCDT a,const char * sch,double sVal,const char * opt)1622 void MGL_EXPORT mgl_contf3_xyz(HMGL gr, HCDT x, HCDT y, HCDT z, HCDT a, const char *sch, double sVal, const char *opt)
1623 {
1624 mreal r = gr->SaveState(opt);
1625 long Num = mgl_isnan(r)?7:long(r+0.5);
1626 if(Num<1) { gr->SetWarn(mglWarnCnt,"ContF3"); return; }
1627 mglDataV v(Num+2); v.Fill(gr->Min.c, gr->Max.c);
1628 mgl_contf3_xyz_val(gr,&v,x,y,z,a,sch,sVal,0);
1629 }
1630 //-----------------------------------------------------------------------------
mgl_contf3(HMGL gr,HCDT a,const char * sch,double sVal,const char * opt)1631 void MGL_EXPORT mgl_contf3(HMGL gr, HCDT a, const char *sch, double sVal, const char *opt)
1632 {
1633 mreal r = gr->SaveState(opt);
1634 long Num = mgl_isnan(r)?7:long(r+0.5);
1635 if(Num<1) { gr->SetWarn(mglWarnCnt,"ContF3"); return; }
1636 mglDataV v(Num+2); v.Fill(gr->Min.c, gr->Max.c);
1637 mgl_contf3_val(gr,&v,a,sch,sVal,0);
1638 }
1639 //-----------------------------------------------------------------------------
mgl_contf3_xyz_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1640 void MGL_EXPORT mgl_contf3_xyz_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1641 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1642 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1643 mgl_contf3_xyz_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1644 delete []o; delete []s; }
1645 //-----------------------------------------------------------------------------
mgl_contf3_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1646 void MGL_EXPORT mgl_contf3_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1647 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1648 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1649 mgl_contf3_val(_GR_, _DA_(v), _DA_(a), s, *sVal, o);
1650 delete []o; delete []s; }
1651 //-----------------------------------------------------------------------------
mgl_contf3_xyz_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * z,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1652 void MGL_EXPORT mgl_contf3_xyz_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *z, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1653 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1654 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1655 mgl_contf3_xyz(_GR_, _DA_(x), _DA_(y), _DA_(z), _DA_(a), s, *sVal, o);
1656 delete []o; delete []s; }
1657 //-----------------------------------------------------------------------------
mgl_contf3_(uintptr_t * gr,uintptr_t * a,const char * sch,mreal * sVal,const char * opt,int l,int lo)1658 void MGL_EXPORT mgl_contf3_(uintptr_t *gr, uintptr_t *a, const char *sch, mreal *sVal, const char *opt,int l,int lo)
1659 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1660 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1661 mgl_contf3(_GR_, _DA_(a), s, *sVal, o);
1662 delete []o; delete []s; }
1663 //-----------------------------------------------------------------------------
1664 //
1665 // Axial series
1666 //
1667 //-----------------------------------------------------------------------------
mgl_find_prev(long i,long pc,long * nn)1668 long MGL_LOCAL_PURE mgl_find_prev(long i, long pc, long *nn)
1669 {
1670 for(long k=0;k<pc;k++) if(nn[k]==i) return k;
1671 return -1;
1672 }
mgl_axial_plot(mglBase * gr,long pc,mglPoint * ff,long * nn,char dir,mreal cc,int wire)1673 void static mgl_axial_plot(mglBase *gr,long pc, mglPoint *ff, long *nn,char dir,mreal cc,int wire)
1674 {
1675 mglPoint a(0,0,1),b,c,q1,q2;
1676 if(dir=='x') a.Set(1,0,0);
1677 if(dir=='y') a.Set(0,1,0);
1678 b = !a; c = a^b;
1679
1680 gr->Reserve(pc*82);
1681 for(long i=0;i<pc;i++)
1682 {
1683 if(nn[i]<0) continue;
1684 long k = mgl_find_prev(i,pc,nn);
1685 q1 = k<0 ? ff[nn[i]]-ff[i] : (ff[nn[i]]-ff[k])*0.5;
1686 q2 = nn[nn[i]]<0 ? ff[nn[i]]-ff[i] : (ff[nn[nn[i]]]-ff[i])*0.5;
1687
1688 long kq = gr->AllocPnts(41*2);
1689 #pragma omp parallel for
1690 for(long j=0;j<41;j++)
1691 {
1692 float co = mgl_cos[(j*18)%360], si = mgl_cos[(270+j*18)%360];
1693 // fi = j*M_PI/20; si = sin(fi); co = cos(fi);
1694 mglPoint p1 = a*ff[i].y + b*(si*ff[i].x) + c*(co*ff[i].x);
1695 mglPoint p2 = a*ff[nn[i]].y + b*(si*ff[nn[i]].x) + c*(co*ff[nn[i]].x);
1696 if(wire)
1697 { gr->AddPntQ(kq+2*j,p1,cc); gr->AddPntQ(kq+2*j+1,p2,cc); }
1698 else
1699 {
1700 gr->AddPntQ(kq+2*j, p1,cc,(a*q1.y + b*(si*q1.x) + c*(co*q1.x))^(b*co-c*si));
1701 gr->AddPntQ(kq+2*j+1,p2,cc,(a*q2.y + b*(si*q2.x) + c*(co*q2.x))^(b*co-c*si));
1702 }
1703 }
1704 if(wire==1)
1705 {
1706 gr->line_plot(kq,kq+1);
1707 for(long j=1;j<41;j++)
1708 {
1709 gr->line_plot(kq+2*j,kq+2*j+1); gr->line_plot(kq+2*j,kq+2*j-2);
1710 gr->line_plot(kq+2*j-1,kq+2*j+1); gr->line_plot(kq+2*j-1,kq+2*j-2);
1711 }
1712 }
1713 else if(wire) for(long j=0;j<41;j++)
1714 { gr->mark_plot(kq+2*j,'.'); gr->mark_plot(kq+2*j+1,'.'); }
1715 else for(long j=1;j<41;j++)
1716 { long i0 = kq+2*j; gr->quad_plot(i0-2,i0-1,i0,i0+1); }
1717 }
1718 }
1719 //-----------------------------------------------------------------------------
1720 // NOTE! All data MUST have the same size! Only first slice is used!
mgl_axial_genI(HMGL gr,mreal val,HCDT a,HCDT x,HCDT y,mreal c,char dir,long ak,int wire)1721 void MGL_EXPORT mgl_axial_genI(HMGL gr, mreal val, HCDT a, HCDT x, HCDT y, mreal c, char dir,long ak,int wire)
1722 {
1723 long n=a->GetNx(), m=a->GetNy();
1724 if(n<2 || m<2 || x->GetNx()*x->GetNy()!=n*m || y->GetNx()*y->GetNy()!=n*m)
1725 { gr->SetWarn(mglWarnDim,"ContGen"); return; }
1726
1727 mglPoint *kk = new mglPoint[2*n*m],*pp = new mglPoint[2*n*m],p;
1728 long pc=0;
1729 // Usually number of points is much smaller. So, there is no reservation.
1730 // gr->Reserve(2*n*m);
1731
1732 // add intersection point of isoline and X or Y axis
1733 const mglData *mx = dynamic_cast<const mglData *>(x);
1734 const mglData *my = dynamic_cast<const mglData *>(y);
1735 const mglData *ma = dynamic_cast<const mglData *>(a);
1736 if(mx&&my&&ma) for(long j=0;j<m;j++) for(long i=0;i<n;i++)
1737 {
1738 long i0 = i+n*j;
1739 mreal d = (i<n-1)?mgl_d(val,ma->a[i0+n*m*ak],ma->a[i0+1+n*m*ak]):-1;
1740 if(d>=0 && d<1)
1741 {
1742 pp[pc].Set(mx->a[i0]*(1-d)+mx->a[i0+1]*d, my->a[i0]*(1-d)+my->a[i0+1]*d);
1743 kk[pc].Set(i+d,j); pc++;
1744 }
1745 d = (j<m-1)?mgl_d(val,ma->a[i0+n*m*ak],ma->a[i0+n*m*ak+n]):-1;
1746 if(d>=0 && d<1)
1747 {
1748 pp[pc].Set(mx->a[i0]*(1-d)+mx->a[i0+n]*d, my->a[i0]*(1-d)+my->a[i0+n]*d);
1749 kk[pc].Set(i,j+d); pc++;
1750 }
1751 }
1752 else for(long j=0;j<m;j++) for(long i=0;i<n;i++)
1753 {
1754 mreal va=a->v(i,j,ak),vx=x->v(i,j),vy=y->v(i,j);
1755 mreal d = (i<n-1)?mgl_d(val,va,a->v(i+1,j,ak)):-1;
1756 if(d>=0 && d<1)
1757 {
1758 pp[pc].Set(vx*(1-d)+x->v(i+1,j)*d, vy*(1-d)+y->v(i+1,j)*d);
1759 kk[pc].Set(i+d,j); pc++;
1760 }
1761 d = (j<m-1)?mgl_d(val,va,a->v(i,j+1,ak)):-1;
1762 if(d>=0 && d<1)
1763 {
1764 pp[pc].Set(vx*(1-d)+x->v(i,j+1)*d, vy*(1-d)+y->v(i,j+1)*d);
1765 kk[pc].Set(i,j+d); pc++;
1766 }
1767 }
1768 // deallocate arrays and finish if no point
1769 if(pc==0) { delete []kk; delete []pp; return; }
1770 // allocate arrays for curve
1771 long *nn = new long[pc], *ff = new long[pc];
1772 for(long i=0;i<pc;i++) nn[i] = ff[i] = -1;
1773 // connect points to line
1774 long j=-1; // current point
1775 do{
1776 if(j>=0)
1777 {
1778 mreal kx = kk[j].x, ky = kk[j].y; long i = -1;
1779 long i11 = long(kx+1e-5), i12 = long(kx-1e-5);
1780 long j11 = long(ky+1e-5), j12 = long(ky-1e-5);
1781 for(long k=0;k<pc;k++) // find closest point in grid
1782 {
1783 if(k==j || k==ff[j] || ff[k]!=-1) continue; // point is marked
1784 long i21 = long(kk[k].x+1e-5), i22 = long(kk[k].x-1e-5);
1785 long j21 = long(kk[k].y+1e-5), j22 = long(kk[k].y-1e-5);
1786 // check if in the same cell
1787 bool cond = (i11==i21 || i11==i22 || i12==i21 || i12==i22) &&
1788 (j11==j21 || j11==j22 || j12==j21 || j12==j22);
1789 if(cond){ i=k; break; }
1790 }
1791 if(i<0) j = -1; // no free close points
1792 else // mark the point
1793 { nn[j] = i; ff[i] = j; j = nn[i]<0 ? i : -1; }
1794 }
1795 if(j<0)
1796 {
1797 for(long k=0;k<pc;k++) if(nn[k]==-1) // first check edges
1798 {
1799 if(kk[k].x==0 || fabs(kk[k].x-n+1)<1e-5 || kk[k].y==0 || fabs(kk[k].y-m+1)<1e-5)
1800 { nn[k]=-2; j = k; break; }
1801 }
1802 if(j<0) for(long k=0;k<pc;k++) if(nn[k]==-1) // or any points inside
1803 { j = k; nn[k]=-2; break; }
1804 }
1805 }while(j>=0);
1806 mgl_axial_plot(gr,pc,pp,nn,dir,c,wire);
1807 delete []kk; delete []nn; delete []ff; delete []pp;
1808 }
1809 //-----------------------------------------------------------------------------
mgl_axial_xy_val(HMGL gr,HCDT v,HCDT x,HCDT y,HCDT z,const char * sch,const char * opt)1810 void MGL_EXPORT mgl_axial_xy_val(HMGL gr, HCDT v, HCDT x, HCDT y, HCDT z, const char *sch, const char *opt)
1811 {
1812 long n=z->GetNx(),m=z->GetNy();
1813 if(mgl_check_dim2(gr,x,y,z,0,"Axial")) return;
1814 gr->SaveState(opt);
1815 static int cgid=1; gr->StartGroup("Axial",cgid++);
1816 long s=gr->AddTexture(sch);
1817 char dir='y';
1818 if(mglchr(sch,'x')) dir = 'x';
1819 if(mglchr(sch,'z')) dir = 'z';
1820
1821 mglData xx, yy;
1822 if(x->GetNx()*x->GetNy()!=m*n || y->GetNx()*y->GetNy()!=m*n) // make
1823 {
1824 xx.Create(n, m); yy.Create(n, m);
1825 for(long i=0;i<n;i++) xx.a[i]=x->v(i);
1826 for(long j=1;j<m;j++) memcpy(xx.a+n*j,xx.a,n*sizeof(mreal));
1827 for(long j=0;j<m;j++)
1828 { mreal t=y->v(j); for(long i=0;i<n;i++) yy.a[i+n*j]=t; }
1829 x = &xx; y = &yy;
1830 }
1831 // x, y -- have the same size z
1832 int wire = mglchr(sch,'#')?1:0;
1833 if(mglchr(sch,'.')) wire = 2;
1834 #pragma omp parallel for collapse(2)
1835 for(long i=0;i<v->GetNx();i++) for(long j=0;j<z->GetNz();j++)
1836 {
1837 if(gr->NeedStop()) continue;
1838 mreal v0 = v->v(i);
1839 mgl_axial_genI(gr,v0,z,x,y,gr->GetC(s,v0),dir,j,wire);
1840 }
1841 gr->EndGroup();
1842 }
1843 //-----------------------------------------------------------------------------
mgl_axial_val(HMGL gr,HCDT v,HCDT a,const char * sch,const char * opt)1844 void MGL_EXPORT mgl_axial_val(HMGL gr, HCDT v, HCDT a, const char *sch, const char *opt)
1845 {
1846 long n=a->GetNx(), m=a->GetNy();
1847 if(n<2 || m<2) { gr->SetWarn(mglWarnLow,"Axial"); return; }
1848 gr->SaveState(opt);
1849 mglDataV x(n, m), y(n, m);
1850 if(gr->Max.x*gr->Min.x>=0) x.Fill(gr->Min.x,gr->Max.x,'x');
1851 else x.Fill(0,gr->Max.x,'x');
1852 y.Fill(gr->Min.y,gr->Max.y,'y');
1853 mgl_axial_xy_val(gr,v,&x,&y,a,sch,0);
1854 }
1855 //-----------------------------------------------------------------------------
mgl_axial_xy(HMGL gr,HCDT x,HCDT y,HCDT a,const char * sch,const char * opt)1856 void MGL_EXPORT mgl_axial_xy(HMGL gr, HCDT x, HCDT y, HCDT a, const char *sch, const char *opt)
1857 {
1858 mreal r = gr->SaveState(opt);
1859 long Num = mgl_isnan(r)?3:long(r+0.5);
1860 if(Num<1) { gr->SetWarn(mglWarnCnt,"Axial"); return; }
1861 mglData v(Num);
1862 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1863 mgl_axial_xy_val(gr,&v,x,y,a,sch,0);
1864 }
1865 //-----------------------------------------------------------------------------
mgl_axial(HMGL gr,HCDT a,const char * sch,const char * opt)1866 void MGL_EXPORT mgl_axial(HMGL gr, HCDT a, const char *sch, const char *opt)
1867 {
1868 mreal r = gr->SaveState(opt);
1869 long Num = mgl_isnan(r)?3:long(r+0.5);
1870 if(Num<1) { gr->SetWarn(mglWarnCnt,"Axial"); return; }
1871 mglData v(Num);
1872 for(long i=0;i<Num;i++) v.a[i] = gr->Min.c + (gr->Max.c-gr->Min.c)*mreal(i+1)/(Num+1);
1873 mgl_axial_val(gr,&v,a,sch,0);
1874 }
1875 //-----------------------------------------------------------------------------
mgl_axial_xy_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1876 void MGL_EXPORT mgl_axial_xy_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1877 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1878 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1879 mgl_axial_xy_val(_GR_, _DA_(v), _DA_(x), _DA_(y), _DA_(a), s, o);
1880 delete []o; delete []s; }
1881 //-----------------------------------------------------------------------------
mgl_axial_val_(uintptr_t * gr,uintptr_t * v,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1882 void MGL_EXPORT mgl_axial_val_(uintptr_t *gr, uintptr_t *v, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1883 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1884 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1885 mgl_axial_val(_GR_, _DA_(v), _DA_(a), s, o); delete []o; delete []s; }
1886 //-----------------------------------------------------------------------------
mgl_axial_xy_(uintptr_t * gr,uintptr_t * x,uintptr_t * y,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1887 void MGL_EXPORT mgl_axial_xy_(uintptr_t *gr, uintptr_t *x, uintptr_t *y, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1888 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1889 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1890 mgl_axial_xy(_GR_, _DA_(x), _DA_(y), _DA_(a), s, o); delete []o; delete []s; }
1891 //-----------------------------------------------------------------------------
mgl_axial_(uintptr_t * gr,uintptr_t * a,const char * sch,const char * opt,int l,int lo)1892 void MGL_EXPORT mgl_axial_(uintptr_t *gr, uintptr_t *a, const char *sch, const char *opt,int l,int lo)
1893 { char *s=new char[l+1]; memcpy(s,sch,l); s[l]=0;
1894 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1895 mgl_axial(_GR_, _DA_(a), s, o); delete []o; delete []s; }
1896 //-----------------------------------------------------------------------------
1897 //
1898 // Torus series
1899 //
1900 //-----------------------------------------------------------------------------
mgl_torus(HMGL gr,HCDT r,HCDT z,const char * sch,const char * opt)1901 void MGL_EXPORT mgl_torus(HMGL gr, HCDT r, HCDT z, const char *sch, const char *opt)
1902 {
1903 long i,j,n=r->GetNx();
1904 if(n*r->GetNy()!=z->GetNx()*z->GetNy()) { gr->SetWarn(mglWarnDim,"Torus"); return; }
1905 if(n<2) { gr->SetWarn(mglWarnLow,"Torus"); return; }
1906 gr->SaveState(opt);
1907 static int cgid=1; gr->StartGroup("Torus",cgid++);
1908
1909 mglPoint *pp = new mglPoint[n];
1910 long *nn = new long[n];
1911 long ss=gr->AddTexture(sch);
1912 char dir='y';
1913 if(mglchr(sch,'x')) dir = 'x';
1914 if(mglchr(sch,'z')) dir = 'z';
1915
1916 mreal c = gr->GetC(ss,gr->Min.c);
1917 const mglData *mr = dynamic_cast<const mglData *>(r);
1918 const mglData *mz = dynamic_cast<const mglData *>(z);
1919 int wire = mglchr(sch,'#')?1:0;
1920 if(mglchr(sch,'.')) wire = 2;
1921 for(j=0;j<r->GetNy();j++)
1922 {
1923 if(mr&&mz) for(i=0;i<n;i++)
1924 {
1925 nn[i] = i<n-1 ? i+1 : -1;
1926 pp[i].Set(mr->a[i+n*j], mz->a[i+n*j]);
1927 }
1928 else for(i=0;i<n;i++)
1929 {
1930 nn[i] = i<n-1 ? i+1 : -1;
1931 pp[i].Set(r->v(i,j), z->v(i,j));
1932 }
1933 mgl_axial_plot(gr,n,pp,nn,dir,c,wire);
1934 }
1935 gr->EndGroup();
1936 delete []nn; delete []pp;
1937 }
1938 //-----------------------------------------------------------------------------
mgl_torus_(uintptr_t * gr,uintptr_t * r,uintptr_t * z,const char * pen,const char * opt,int l,int lo)1939 void MGL_EXPORT mgl_torus_(uintptr_t *gr, uintptr_t *r, uintptr_t *z, const char *pen, const char *opt,int l,int lo)
1940 { char *s=new char[l+1]; memcpy(s,pen,l); s[l]=0;
1941 char *o=new char[lo+1]; memcpy(o,opt,lo); o[lo]=0;
1942 mgl_torus(_GR_, _DA_(r), _DA_(z), s, o); delete []o; delete []s; }
1943 //-----------------------------------------------------------------------------
1944