diff options
Diffstat (limited to 'Master/texmf/asymptote/graph.asy')
-rw-r--r-- | Master/texmf/asymptote/graph.asy | 2214 |
1 files changed, 2214 insertions, 0 deletions
diff --git a/Master/texmf/asymptote/graph.asy b/Master/texmf/asymptote/graph.asy new file mode 100644 index 00000000000..805cf15f612 --- /dev/null +++ b/Master/texmf/asymptote/graph.asy @@ -0,0 +1,2214 @@ +private import math; + +import splinetype; +import graph_settings; + +scaleT Linear; + +scaleT Log=scaleT(log10,pow10,logarithmic=true); +scaleT Logarithmic=Log; + +string baselinetemplate="$10^4$"; + +// A linear scale, with optional autoscaling of minimum and maximum values, +// scaling factor s and intercept. +scaleT Linear(bool automin=false, bool automax=automin, real s=1, + real intercept=0) +{ + real sinv=1/s; + scalefcn T,Tinv; + if(s == 1 && intercept == 0) + T=Tinv=identity; + else { + T=new real(real x) {return (x-intercept)*s;}; + Tinv=new real(real x) {return x*sinv+intercept;}; + } + return scaleT(T,Tinv,logarithmic=false,automin,automax); +} + +// A logarithmic scale, with optional autoscaling of minimum and maximum +// values. +scaleT Log(bool automin=false, bool automax=automin) +{ + return scaleT(Log.T,Log.Tinv,logarithmic=true,automin,automax); +} + +// A "broken" linear axis omitting the segment [a,b]. +scaleT Broken(real a, real b, bool automin=false, bool automax=automin) +{ + real skip=b-a; + real T(real x) { + if(x <= a) return x; + if(x <= b) return a; + return x-skip; + } + real Tinv(real x) { + if(x <= a) return x; + return x+skip; + } + return scaleT(T,Tinv,logarithmic=false,automin,automax); +} + +// A "broken" logarithmic axis omitting the segment [a,b], where a and b are +// automatically rounded to the nearest integral power of the base. +scaleT BrokenLog(real a, real b, bool automin=false, bool automax=automin) +{ + real A=round(Log.T(a)); + real B=round(Log.T(b)); + a=Log.Tinv(A); + b=Log.Tinv(B); + real skip=B-A; + real T(real x) { + if(x <= a) return Log.T(x); + if(x <= b) return A; + return Log.T(x)-skip; + } + real Tinv(real x) { + real X=Log.Tinv(x); + if(X <= a) return X; + return Log.Tinv(x+skip); + } + return scaleT(T,Tinv,logarithmic=true,automin,automax); +} + +Label Break=Label("$\approx$",UnFill(0.2mm)); + +void scale(picture pic=currentpicture, scaleT x, scaleT y=x, scaleT z=y) +{ + pic.scale.x.scale=x; + pic.scale.y.scale=y; + pic.scale.z.scale=z; + pic.scale.x.automin=x.automin; + pic.scale.y.automin=y.automin; + pic.scale.z.automin=z.automin; + pic.scale.x.automax=x.automax; + pic.scale.y.automax=y.automax; + pic.scale.z.automax=z.automax; +} + +void scale(picture pic=currentpicture, bool xautoscale=false, + bool yautoscale=xautoscale, bool zautoscale=yautoscale) +{ + scale(pic,Linear(xautoscale,xautoscale),Linear(yautoscale,yautoscale), + Linear(zautoscale,zautoscale)); +} + +struct scientific +{ + int sign; + real mantissa; + int exponent; + int ceil() {return sign*ceil(mantissa);} + real scale(real x, real exp) { + static real max=0.1*realMax; + static real limit=-log10(max); + return x*(exp > limit ? 10^-exp : max); + } + real ceil(real x, real exp) {return ceil(sign*scale(abs(x),exp));} + real floor(real x, real exp) {return floor(sign*scale(abs(x),exp));} +} + +// Convert x to scientific notation +scientific scientific(real x) +{ + scientific s; + s.sign=sgn(x); + x=abs(x); + if(x == 0) {s.mantissa=0; s.exponent=-intMax; return s;} + real logx=log10(x); + s.exponent=floor(logx); + s.mantissa=s.scale(x,s.exponent); + return s; +} + +// Autoscale limits and tick divisor. +struct bounds { + real min; + real max; + + // Possible tick intervals: + int[] divisor; + + void operator init(real min, real max, int[] divisor=new int[]) { + this.min=min; + this.max=max; + this.divisor=divisor; + } +} + +// Compute tick divisors. +int[] divisors(int a, int b) +{ + int[] dlist; + int n=b-a; + dlist[0]=1; + if(n == 1) {dlist[1]=10; dlist[2]=100; return dlist;} + if(n == 2) {dlist[1]=2; return dlist;} + int sqrtn=floor(sqrt(n)); + int i=0; + for(int d=2; d <= sqrtn; ++d) + if(n % d == 0 && (a*b >= 0 || b % (n/d) == 0)) dlist[++i]=d; + for(int d=sqrtn; d >= 1; --d) + if(n % d == 0 && (a*b >= 0 || b % d == 0)) dlist[++i]=quotient(n,d); + return dlist; +} + +real upscale(real b, real a) +{ + if(b <= 5) b=5; + else if (b > 10 && a >= 0 && b <= 12) b=12; + else if (b > 10 && (a >= 0 || 15 % -a == 0) && b <= 15) b=15; + else b=ceil(b/10)*10; + return b; +} + +// Compute autoscale limits and tick divisor. +bounds autoscale(real Min, real Max, scaleT scale=Linear) +{ + bounds m; + if(scale.logarithmic) { + m.min=floor(Min); + m.max=ceil(Max); + return m; + } + if(!(finite(Min) && finite(Max))) + abort("autoscale requires finite limits"); + Min=scale.Tinv(Min); + Max=scale.Tinv(Max); + m.min=Min; + m.max=Max; + if(Min > Max) {real temp=Min; Min=Max; Max=temp;} + if(Min == Max) { + if(Min == 0) {m.max=1; return m;} + if(Min > 0) {Min=0; Max *= 2;} + else {Min *= 2; Max=0;} + } + + int sign; + if(Min < 0 && Max <= 0) {real temp=-Min; Min=-Max; Max=temp; sign=-1;} + else sign=1; + scientific sa=scientific(Min); + scientific sb=scientific(Max); + int exp=max(sa.exponent,sb.exponent); + real a=sa.floor(Min,exp); + real b=sb.ceil(Max,exp); + + void zoom() { + --exp; + a=sa.floor(Min,exp); + b=sb.ceil(Max,exp); + } + + if(sb.mantissa <= 1.5) + zoom(); + + while((b-a)*10.0^exp > 10*(Max-Min)) + zoom(); + + real bsave=b; + if(b-a > (a >= 0 ? 8 : 6)) { + b=upscale(b,a); + if(a >= 0) { + if(a <= 5) a=0; else a=floor(a/10)*10; + } else a=-upscale(-a,-1); + } + + // Redo b in case the value of a has changed + if(bsave-a > (a >= 0 ? 8 : 6)) + b=upscale(bsave,a); + + if(sign == -1) {real temp=-a; a=-b; b=temp;} + real Scale=10.0^exp; + m.min=scale.T(a*Scale); + m.max=scale.T(b*Scale); + if(m.min > m.max) {real temp=m.min; m.min=m.max; m.max=temp;} + m.divisor=divisors(round(a),round(b)); + return m; +} + +typedef string ticklabel(real); + +ticklabel Format(string s=defaultformat) +{ + return new string(real x) {return format(s,x);}; +} + +ticklabel OmitFormat(string s=defaultformat ... real[] x) +{ + return new string(real v) { + int i=find(x == v); + return i < 0 ? format(s,v) : ""; + }; +} + +string trailingzero="$%#$"; +string signedtrailingzero="$%+#$"; + +ticklabel DefaultFormat=Format(); +ticklabel NoZeroFormat=OmitFormat(0); + +// Format tick values as integral powers of base; otherwise with DefaultFormat. +ticklabel DefaultLogFormat(int base) { + return new string(real x) { + string exponent=format("%.4f",log(x)/log(base)); + return find(exponent,".") == -1 ? + "$"+(string) base+"^{"+exponent+"}$" : format(x); + }; +} + +// Format all tick values as powers of base. +ticklabel LogFormat(int base) { + return new string(real x) { + return format("$"+(string) base+"^{%g}$",log(x)/log(base)); + }; +} + +ticklabel LogFormat=LogFormat(10); +ticklabel DefaultLogFormat=DefaultLogFormat(10); + +// The default direction specifier. +pair zero(real) {return 0;} + +struct ticklocate { + real a,b; // Tick values at point(g,0), point(g,length(g)). + autoscaleT S; // Autoscaling transformation. + pair dir(real t); // Absolute 2D tick direction. + triple dir3(real t); // Absolute 3D tick direction. + real time(real v); // Returns the time corresponding to the value v. + ticklocate copy() { + ticklocate T=new ticklocate; + T.a=a; + T.b=b; + T.S=S.copy(); + T.dir=dir; + T.dir3=dir3; + T.time=time; + return T; + } +} + +autoscaleT defaultS; + +typedef real valuetime(real); + +valuetime linear(scalefcn S=identity, real Min, real Max) +{ + real factor=Max == Min ? 0.0 : 1.0/(Max-Min); + return new real(real v) {return (S(v)-Min)*factor;}; +} + +ticklocate ticklocate(real a, real b, autoscaleT S=defaultS, + real tickmin=-infinity, real tickmax=infinity, + real time(real)=null, pair dir(real)=zero) +{ + if((valuetime) time == null) time=linear(S.T(),a,b); + ticklocate locate; + locate.a=a; + locate.b=b; + locate.S=S.copy(); + if(finite(tickmin)) locate.S.tickMin=tickmin; + if(finite(tickmax)) locate.S.tickMax=tickmax; + locate.time=time; + locate.dir=dir; + return locate; +} + +private struct locateT { + real t; // tick location time + pair Z; // tick location in frame coordinates + pair pathdir; // path direction in frame coordinates + pair dir; // tick direction in frame coordinates + + void dir(transform T, path g, ticklocate locate, real t) { + pathdir=unit(shiftless(T)*dir(g,t)); + pair Dir=locate.dir(t); + dir=Dir == 0 ? -I*pathdir : unit(Dir); + } + // Locate the desired position of a tick along a path. + void calc(transform T, path g, ticklocate locate, real val) { + t=locate.time(val); + Z=T*point(g,t); + dir(T,g,locate,t); + } +} + +pair ticklabelshift(pair align, pen p=currentpen) +{ + return 0.25*unit(align)*labelmargin(p); +} + +void drawtick(frame f, transform T, path g, path g2, ticklocate locate, + real val, real Size, int sign, pen p, bool extend) +{ + locateT locate1,locate2; + locate1.calc(T,g,locate,val); + if(extend && size(g2) > 0) { + locate2.calc(T,g2,locate,val); + draw(f,locate1.Z--locate2.Z,p); + } else + if(sign == 0) + draw(f,locate1.Z-Size*locate1.dir--locate1.Z+Size*locate1.dir,p); + else + draw(f,locate1.Z--locate1.Z+Size*sign*locate1.dir,p); +} + +real zerotickfuzz=10*epsilon; + +// Label a tick on a frame. +pair labeltick(frame d, transform T, path g, ticklocate locate, real val, + pair side, int sign, real Size, ticklabel ticklabel, + Label F, real norm=0) +{ + locateT locate1; + locate1.calc(T,g,locate,val); + pair align=side*locate1.dir; + pair perp=I*locate1.pathdir; + + // Adjust tick label alignment + pair adjust=unit(align+0.75perp*sgn(dot(align,perp))); + // Project align onto adjusted direction. + align=adjust*dot(align,adjust); + pair shift=dot(align,-sign*locate1.dir) <= 0 ? align*Size : + ticklabelshift(align,F.p); + + real label; + if(locate.S.scale.logarithmic) + label=locate.S.scale.Tinv(val); + else { + label=val; + if(abs(label) < zerotickfuzz*norm) label=0; + // Fix epsilon errors at +/-1e-4 + // default format changes to scientific notation here + if(abs(abs(label)-1e-4) < epsilon) label=sgn(label)*1e-4; + } + + string s=ticklabel(label); + if(s != "") + label(d,F.T*baseline(s,baselinetemplate),locate1.Z+shift,align,F.p, + F.filltype); + return locate1.pathdir; +} + +// Add axis label L to frame f. +void labelaxis(frame f, transform T, Label L, path g, + ticklocate locate=null, int sign=1, bool ticklabels=false) +{ + Label L0=L.copy(); + real t=L0.relative(g); + pair z=point(g,t); + pair dir=dir(g,t); + pair perp=I*dir; + if(locate != null) { + locateT locate1; + locate1.dir(T,g,locate,t); + L0.align(L0.align,unit(-sgn(dot(sign*locate1.dir,perp))*perp)); + } + pair align=L0.align.dir; + if(L0.align.relative) align *= -perp; + pair alignperp=dot(align,perp)*perp; + pair offset; + if(ticklabels) { + if(piecewisestraight(g)) { + real angle=degrees(dir,warn=false); + transform S=rotate(-angle,z); + frame F=S*f; + pair Align=rotate(-angle)*alignperp; + offset=unit(alignperp-sign*locate.dir(t))* + abs((Align.y >= 0 ? max(F).y : (Align.y < 0 ? min(F).y : 0))-z.y); + } + z += offset; + } + + L0.align(align); + L0.position(z); + frame d; + add(d,L0); + pair width=0.5*size(d); + int n=length(g); + real t=L.relative(); + pair s=realmult(width,dir(g,t)); + if(t <= 0) { + if(L.align.default) s *= -axislabelfactor; + d=shift(s)*d; + } else if(t >= n) { + if(L.align.default) s *= -axislabelfactor; + d=shift(-s)*d; + } else if(offset == 0 && L.align.default) { + pair s=realmult(width,I*dir(g,t)); + s=axislabelfactor*s; + d=shift(s)*d; + } + add(f,d); +} + +// Check the tick coverage of a linear axis. +bool axiscoverage(int N, transform T, path g, ticklocate locate, real Step, + pair side, int sign, real Size, Label F, ticklabel ticklabel, + real norm, real limit) +{ + real coverage=0; + bool loop=cyclic(g); + real a=locate.S.Tinv(locate.a); + real b=locate.S.Tinv(locate.b); + real tickmin=finite(locate.S.tickMin) ? locate.S.Tinv(locate.S.tickMin) : a; + if(Size > 0) { + int count=0; + if(loop) count=N+1; + else { + for(int i=0; i <= N; ++i) { + real val=tickmin+i*Step; + if(val >= a && val <= b) + ++count; + } + } + if(count > 0) limit /= count; + for(int i=0; i <= N; ++i) { + real val=tickmin+i*Step; + if(loop || (val >= a && val <= b)) { + frame d; + pair dir=labeltick(d,T,g,locate,val,side,sign,Size,ticklabel,F,norm); + if(abs(dot(size(d),dir)) > limit) return false; + } + } + } + return true; +} + +// Check the tick coverage of a logarithmic axis. +bool logaxiscoverage(int N, transform T, path g, ticklocate locate, pair side, + int sign, real Size, Label F, ticklabel ticklabel, + real limit, int first, int last) +{ + bool loop=cyclic(g); + real coverage=0; + real a=locate.a; + real b=locate.b; + int count=0; + for(int i=first-1; i <= last+1; i += N) { + if(loop || i >= a && i <= b) + ++count; + } + if(count > 0) limit /= count; + for(int i=first-1; i <= last+1; i += N) { + if(loop || i >= a && i <= b) { + frame d; + pair dir=labeltick(d,T,g,locate,i,side,sign,Size,ticklabel,F); + if(abs(dot(size(d),dir)) > limit) return false; + } + } + return true; +} + +struct tickvalues { + real[] major; + real[] minor; + int N; // For logarithmic axes: number of decades between tick labels. +} + +// Determine a format that distinguishes adjacent pairs of ticks, optionally +// adding trailing zeros. +string autoformat(string format="", real norm ... real[] a) +{ + bool trailingzero=(format == trailingzero); + bool signedtrailingzero=(format == signedtrailingzero); + if(!trailingzero && !signedtrailingzero && format != "") return format; + + real[] A=sort(a); + real[] a=abs(A); + + bool signchange=(A.length > 1 && A[0] < 0 && A[A.length-1] >= 0); + + for(int i=0; i < A.length-1; ++i) + if(a[i] < zerotickfuzz*norm) A[i]=a[i]=0; + + int n=0; + + bool Fixed=find(a >= 1e4-epsilon | (a > 0 & a <= 1e-4-epsilon)) < 0; + string Format=defaultformat(4,fixed=Fixed); + + if(Fixed && n < 4) { + for(int i=0; i < A.length-1; ++i) { + real a=A[i]; + while(format(defaultformat(n,fixed=Fixed),a) != format(Format,a)) + ++n; + } + } + + string trailing=trailingzero ? (signchange ? "# " : "#") : + signedtrailingzero ? "#+" : ""; + + string format=defaultformat(n,trailing,Fixed); + + for(int i=0; i < A.length-1; ++i) { + real a=A[i]; + real b=A[i+1]; + // Check if an extra digit of precision should be added. + string fixedformat="%#."+string(n+1)+"f"; + string A=format(fixedformat,a); + string B=format(fixedformat,b); + if(substr(A,length(A)-1,1) != "0" || substr(B,length(B)-1,1) != "0") { + a *= 0.1; + b *= 0.1; + } + if(a != b) { + while(format(format,a) == format(format,b)) + format=defaultformat(++n,trailing,Fixed); + } + } + + if(n == 0) return defaultformat; + return format; +} + +// Automatic tick generation routine. +tickvalues generateticks(int sign, Label F="", ticklabel ticklabel=null, + int N, int n=0, real Step=0, real step=0, + real Size=0, real size=0, + transform T, pair side, path g, real limit, + pen p, ticklocate locate, int[] divisor, + bool opposite) +{ + tickvalues tickvalues; + sign=opposite ? -sign : sign; + if(Size == 0) Size=Ticksize; + if(size == 0) size=ticksize; + F=F.copy(); + F.p(p); + + if(F.align.dir != 0) side=F.align.dir; + else if(side == 0) side=((sign == 1) ? left : right); + + bool ticklabels=false; + path G=T*g; + + if(!locate.S.scale.logarithmic) { + real a=locate.S.Tinv(locate.a); + real b=locate.S.Tinv(locate.b); + real norm=max(abs(a),abs(b)); + string format=autoformat(F.s,norm,a,b); + if(F.s == "%") F.s=""; + if(ticklabel == null) ticklabel=Format(format); + + if(a > b) {real temp=a; a=b; b=temp;} + + if(b-a < 100.0*epsilon*norm) b=a; + + real tickmin=finite(locate.S.tickMin) && (Step == 0 || locate.S.automin) ? + locate.S.Tinv(locate.S.tickMin) : a; + real tickmax=finite(locate.S.tickMax) && (Step == 0 || locate.S.automax) ? + locate.S.Tinv(locate.S.tickMax) : b; + if(tickmin > tickmax) {real temp=tickmin; tickmin=tickmax; tickmax=temp;} + + real inStep=Step; + + bool calcStep=true; + real len=tickmax-tickmin; + if(Step == 0 && N == 0) { + N=1; + if(divisor.length > 0) { + bool autoscale=locate.S.automin && locate.S.automax; + real h=0.5*(b-a); + for(int d=divisor.length-1; d >= 0; --d) { + int N0=divisor[d]; + Step=len/N0; + int N1=N0; + int m=2; + while(Step > h) { + N0=m*N1; + Step=len/N0; + m *= 2; + } + if(axiscoverage(N0,T,g,locate,Step,side,sign,Size,F,ticklabel,norm, + limit)) { + N=N0; + if(N0 == 1 && !autoscale && d < divisor.length-1) { + // Try using 2 ticks (otherwise 1); + int div=divisor[d+1]; + Step=quotient(div,2)*len/div; + calcStep=false; + if(axiscoverage(2,T,g,locate,Step,side,sign,Size,F,ticklabel, + norm,limit)) N=2; + else Step=len; + } + // Found a good divisor; now compute subtick divisor + if(n == 0) { + if(step != 0) n=ceil(Step/step); + else { + n=quotient(divisor[divisor.length-1],N); + if(N == 1) n=(a*b >= 0) ? 2 : 1; + if(n == 1) n=2; + } + } + break; + } + } + } + } + + if(inStep != 0 && !locate.S.automin) { + tickmin=floor(tickmin/Step)*Step; + len=tickmax-tickmin; + } + + if(calcStep) { + if(N == 0) N=(int) (len/Step); + else Step=len/N; + } + + if(n == 0) { + if(step != 0) n=ceil(Step/step); + } else step=Step/n; + + b += epsilon*norm; + + if(Size > 0) { + for(int i=0; i <= N; ++i) { + real val=tickmin+i*Step; + if(val >= a && val <= b) + tickvalues.major.push(val); + if(size > 0 && step > 0) { + real iStep=i*Step; + real jstop=(len-iStep)/step; + for(int j=1; j < n && j <= jstop; ++j) { + real val=tickmin+iStep+j*step; + if(val >= a && val <= b) + tickvalues.minor.push(val); + } + } + } + } + + } else { // Logarithmic + string format=F.s; + if(F.s == "%") F.s=""; + + int base=round(locate.S.scale.Tinv(1)); + + if(ticklabel == null) + ticklabel=format == "%" ? Format("") : DefaultLogFormat(base); + real a=locate.S.postscale.Tinv(locate.a); + real b=locate.S.postscale.Tinv(locate.b); + if(a > b) {real temp=a; a=b; b=temp;} + + int first=floor(a-epsilon); + int last=ceil(b+epsilon); + + if(N == 0) { + N=1; + while(N <= last-first) { + if(logaxiscoverage(N,T,g,locate,side,sign,Size,F,ticklabel,limit, + first,last)) break; + ++N; + } + } + + if(N <= 2 && n == 0) n=base; + tickvalues.N=N; + + if(N > 0) { + for(int i=first-1; i <= last+1; ++i) { + if(i >= a && i <= b) + tickvalues.major.push(locate.S.scale.Tinv(i)); + if(n > 0) { + for(int j=2; j < n; ++j) { + real val=(i+1+locate.S.scale.T(j/n)); + if(val >= a && val <= b) + tickvalues.minor.push(locate.S.scale.Tinv(val)); + } + } + } + } + } + return tickvalues; +} + +// Signature of routines that draw labelled paths with ticks and tick labels. +typedef void ticks(frame, transform, Label, pair, path, path, pen, + arrowbar, margin, ticklocate, int[], bool opposite=false); + +// Tick construction routine for a user-specified array of tick values. +ticks Ticks(int sign, Label F="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + real[] Ticks=new real[], real[] ticks=new real[], int N=1, + bool begin=true, bool end=true, + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return new void(frame f, transform t, Label L, pair side, path g, path g2, + pen p, arrowbar arrow, margin margin, ticklocate locate, + int[] divisor, bool opposite) { + // Use local copy of context variables: + int sign=opposite ? -sign : sign; + pen pTick=pTick; + pen ptick=ptick; + ticklabel ticklabel=ticklabel; + + real Size=Size; + real size=size; + if(Size == 0) Size=Ticksize; + if(size == 0) size=ticksize; + + Label L=L.copy(); + Label F=F.copy(); + L.p(p); + F.p(p); + if(pTick == nullpen) pTick=p; + if(ptick == nullpen) ptick=pTick; + + if(F.align.dir != 0) side=F.align.dir; + else if(side == 0) side=F.T*((sign == 1) ? left : right); + + bool ticklabels=false; + path G=t*g; + path G2=t*g2; + + scalefcn T; + + real a,b; + if(locate.S.scale.logarithmic) { + a=locate.S.postscale.Tinv(locate.a); + b=locate.S.postscale.Tinv(locate.b); + T=locate.S.scale.T; + } else { + a=locate.S.Tinv(locate.a); + b=locate.S.Tinv(locate.b); + T=identity; + } + + if(a > b) {real temp=a; a=b; b=temp;} + + real norm=max(abs(a),abs(b)); + + string format=autoformat(F.s,norm...Ticks); + if(F.s == "%") F.s=""; + if(ticklabel == null) { + if(locate.S.scale.logarithmic) { + int base=round(locate.S.scale.Tinv(1)); + ticklabel=format == "%" ? Format("") : DefaultLogFormat(base); + } else ticklabel=Format(format); + } + + begingroup(f); + if(opposite) draw(f,G,p); + else draw(f,margin(G,p).g,p,arrow); + for(int i=(begin ? 0 : 1); i < (end ? Ticks.length : Ticks.length-1); ++i) { + real val=T(Ticks[i]); + if(val >= a && val <= b) + drawtick(f,t,g,g2,locate,val,Size,sign,pTick,extend); + } + for(int i=0; i < ticks.length; ++i) { + real val=T(ticks[i]); + if(val >= a && val <= b) + drawtick(f,t,g,g2,locate,val,size,sign,ptick,extend); + } + endgroup(f); + + if(N == 0) N=1; + if(Size > 0 && !opposite) { + for(int i=(beginlabel ? 0 : 1); + i < (endlabel ? Ticks.length : Ticks.length-1); i += N) { + real val=T(Ticks[i]); + if(val >= a && val <= b) { + ticklabels=true; + labeltick(f,t,g,locate,val,side,sign,Size,ticklabel,F,norm); + } + } + } + if(L.s != "" && !opposite) + labelaxis(f,t,L,G,locate,sign,ticklabels); + }; +} + +// Optional routine to allow modification of auto-generated tick values. +typedef tickvalues tickmodifier(tickvalues); +tickvalues None(tickvalues v) {return v;} + +tickmodifier OmitTick(... real[] x) { + return new tickvalues(tickvalues v) { + if(v.major.length == 0) return v; + real norm=max(abs(v.major)); + for(int i=0; i < x.length; ++i) { + int j=find(abs(v.major-x[i]) < zerotickfuzz*norm); + if(j >= 0) v.major.delete(j); + } + return v; + }; +} + +tickmodifier NoZero=OmitTick(0); + +// Tickmodifier that removes all major ticks in the interval [a,b]. +tickmodifier Break(real a, real b) { + return new tickvalues(tickvalues v) { + real[] V=v.major; + real[] major; + for(int i=0; i < V.length; ++i) + if(V[i] < a-epsilon || V[i] > b+epsilon) major.push(V[i]); + v.major=major; + return v; + }; +} + +// Automatic tick construction routine. +ticks Ticks(int sign, Label F="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + int N, int n=0, real Step=0, real step=0, + bool begin=true, bool end=true, tickmodifier modify=None, + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return new void(frame f, transform T, Label L, pair side, path g, path g2, + pen p, arrowbar arrow, margin margin, ticklocate locate, + int[] divisor, bool opposite) { + real limit=Step == 0 ? axiscoverage*arclength(T*g) : 0; + tickvalues values=modify(generateticks(sign,F,ticklabel,N,n,Step,step, + Size,size,T,side,g, + limit,p,locate,divisor,opposite)); + + Ticks(sign,F,ticklabel,beginlabel,endlabel,values.major,values.minor, + values.N,begin,end,Size,size,extend,pTick,ptick) + (f,T,L,side,g,g2,p,arrow,margin,locate,divisor,opposite); + }; +} + +ticks NoTicks() +{ + return new void(frame f, transform T, Label L, pair, path g, path, pen p, + arrowbar arrow, margin margin, ticklocate, + int[], bool opposite) { + path G=T*g; + if(opposite) draw(f,G,p); + else { + draw(f,margin(G,p).g,p,arrow); + if(L.s != "") { + Label L=L.copy(); + L.p(p); + labelaxis(f,T,L,G); + } + } + }; +} + +ticks LeftTicks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + int N=0, int n=0, real Step=0, real step=0, + bool begin=true, bool end=true, tickmodifier modify=None, + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(-1,format,ticklabel,beginlabel,endlabel,N,n,Step,step, + begin,end,modify,Size,size,extend,pTick,ptick); +} + +ticks RightTicks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + int N=0, int n=0, real Step=0, real step=0, + bool begin=true, bool end=true, tickmodifier modify=None, + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(1,format,ticklabel,beginlabel,endlabel,N,n,Step,step, + begin,end,modify,Size,size,extend,pTick,ptick); +} + +ticks Ticks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + int N=0, int n=0, real Step=0, real step=0, + bool begin=true, bool end=true, tickmodifier modify=None, + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(0,format,ticklabel,beginlabel,endlabel,N,n,Step,step, + begin,end,modify,Size,size,extend,pTick,ptick); +} + +ticks LeftTicks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + real[] Ticks, real[] ticks=new real[], + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(-1,format,ticklabel,beginlabel,endlabel, + Ticks,ticks,Size,size,extend,pTick,ptick); +} + +ticks RightTicks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + real[] Ticks, real[] ticks=new real[], + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(1,format,ticklabel,beginlabel,endlabel, + Ticks,ticks,Size,size,extend,pTick,ptick); +} + +ticks Ticks(Label format="", ticklabel ticklabel=null, + bool beginlabel=true, bool endlabel=true, + real[] Ticks, real[] ticks=new real[], + real Size=0, real size=0, bool extend=false, + pen pTick=nullpen, pen ptick=nullpen) +{ + return Ticks(0,format,ticklabel,beginlabel,endlabel, + Ticks,ticks,Size,size,extend,pTick,ptick); +} + +ticks NoTicks=NoTicks(), +LeftTicks=LeftTicks(), +RightTicks=RightTicks(), +Ticks=Ticks(); + +pair tickMin(picture pic) +{ + return minbound(pic.userMin(),(pic.scale.x.tickMin,pic.scale.y.tickMin)); +} + +pair tickMax(picture pic) +{ + return maxbound(pic.userMax(),(pic.scale.x.tickMax,pic.scale.y.tickMax)); +} + +int Min=-1; +int Value=0; +int Max=1; +int Both=2; + +// Structure used to communicate axis and autoscale settings to tick routines. +struct axisT { + int type; // -1 = min, 0 = given value, 1 = max, 2 = min/max + int type2; // for 3D axis + real value; + real value2; + pair side; // 2D tick label direction relative to path (left or right) + real position; // label position along axis + align align; // default axis label alignment and 3D tick label direction + int[] xdivisor; + int[] ydivisor; + int[] zdivisor; + bool extend; // extend axis to graph boundary? +}; + +axisT axis; +typedef void axis(picture, axisT); + +axis Bottom(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Min; + axis.position=0.5; + axis.side=right; + axis.align=S; + axis.extend=extend; + }; +} + +axis Top(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Max; + axis.position=0.5; + axis.side=left; + axis.align=N; + axis.extend=extend; + }; +} + +axis BottomTop(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Both; + axis.position=0.5; + axis.side=right; + axis.align=S; + axis.extend=extend; + }; +} + +axis Left(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Min; + axis.position=0.5; + axis.side=left; + axis.align=W; + axis.extend=extend; + }; +} + +axis Right(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Max; + axis.position=0.5; + axis.side=right; + axis.align=E; + axis.extend=extend; + }; +} + +axis LeftRight(bool extend=false) +{ + return new void(picture pic, axisT axis) { + axis.type=Both; + axis.position=0.5; + axis.side=left; + axis.align=W; + axis.extend=extend; + }; +} + +axis XEquals(real x, bool extend=true) +{ + return new void(picture pic, axisT axis) { + axis.type=Value; + axis.value=pic.scale.x.T(x); + axis.position=1; + axis.side=left; + axis.align=W; + axis.extend=extend; + }; +} + +axis YEquals(real y, bool extend=true) +{ + return new void(picture pic, axisT axis) { + axis.type=Value; + axis.value=pic.scale.y.T(y); + axis.position=1; + axis.side=right; + axis.align=S; + axis.extend=extend; + }; +} + +axis XZero(bool extend=true) +{ + return new void(picture pic, axisT axis) { + axis.type=Value; + axis.value=pic.scale.x.T(pic.scale.x.scale.logarithmic ? 1 : 0); + axis.position=1; + axis.side=left; + axis.align=W; + axis.extend=extend; + }; +} + +axis YZero(bool extend=true) +{ + return new void(picture pic, axisT axis) { + axis.type=Value; + axis.value=pic.scale.y.T(pic.scale.y.scale.logarithmic ? 1 : 0); + axis.position=1; + axis.side=right; + axis.align=S; + axis.extend=extend; + }; +} + +axis Bottom=Bottom(), +Top=Top(), +BottomTop=BottomTop(), +Left=Left(), +Right=Right(), +LeftRight=LeftRight(), +XZero=XZero(), +YZero=YZero(); + +// Draw a general axis. +void axis(picture pic=currentpicture, Label L="", path g, path g2=nullpath, + pen p=currentpen, ticks ticks, ticklocate locate, + arrowbar arrow=None, margin margin=NoMargin, + int[] divisor=new int[], bool above=false, bool opposite=false) +{ + Label L=L.copy(); + real t=reltime(g,0.5); + if(L.defaultposition) L.position(t); + divisor=copy(divisor); + locate=locate.copy(); + pic.add(new void (frame f, transform t, transform T, pair lb, pair rt) { + frame d; + ticks(d,t,L,0,g,g2,p,arrow,margin,locate,divisor,opposite); + (above ? add : prepend)(f,t*T*inverse(t)*d); + }); + + pic.addPath(g,p); + + if(L.s != "") { + frame f; + Label L0=L.copy(); + L0.position(0); + add(f,L0); + pair pos=point(g,L.relative()*length(g)); + pic.addBox(pos,pos,min(f),max(f)); + } +} + +real xtrans(transform t, real x) +{ + return (t*(x,0)).x; +} + +real ytrans(transform t, real y) +{ + return (t*(0,y)).y; +} + +// An internal routine to draw an x axis at a particular y value. +void xaxisAt(picture pic=currentpicture, Label L="", axis axis, + real xmin=-infinity, real xmax=infinity, pen p=currentpen, + ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin, + bool above=true, bool opposite=false) +{ + real y=axis.value; + real y2; + Label L=L.copy(); + int[] divisor=copy(axis.xdivisor); + pair side=axis.side; + int type=axis.type; + + pic.add(new void (frame f, transform t, transform T, pair lb, pair rt) { + transform tinv=inverse(t); + pair a=xmin == -infinity ? tinv*(lb.x-min(p).x,ytrans(t,y)) : (xmin,y); + pair b=xmax == infinity ? tinv*(rt.x-max(p).x,ytrans(t,y)) : (xmax,y); + pair a2=xmin == -infinity ? tinv*(lb.x-min(p).x,ytrans(t,y2)) : (xmin,y2); + pair b2=xmax == infinity ? tinv*(rt.x-max(p).x,ytrans(t,y2)) : (xmax,y2); + + if(xmin == -infinity || xmax == infinity) { + bounds mx=autoscale(a.x,b.x,pic.scale.x.scale); + pic.scale.x.tickMin=mx.min; + pic.scale.x.tickMax=mx.max; + divisor=mx.divisor; + } + + real fuzz=epsilon*max(abs(a.x),abs(b.x)); + a -= (fuzz,0); + b += (fuzz,0); + + frame d; + ticks(d,t,L,side,a--b,finite(y2) ? a2--b2 : nullpath,p,arrow,margin, + ticklocate(a.x,b.x,pic.scale.x),divisor,opposite); + (above ? add : prepend)(f,t*T*tinv*d); + }); + + void bounds() { + if(type == Min) + y=pic.scale.y.automin() ? tickMin(pic).y : pic.userMin.y; + else if(type == Max) + y=pic.scale.y.automax() ? tickMax(pic).y : pic.userMax.y; + else if(type == Both) { + y2=pic.scale.y.automax() ? tickMax(pic).y : pic.userMax.y; + y=opposite ? y2 : + (pic.scale.y.automin() ? tickMin(pic).y : pic.userMin.y); + } + + real Xmin=finite(xmin) ? xmin : pic.userMin.x; + real Xmax=finite(xmax) ? xmax : pic.userMax.x; + + pair a=(Xmin,y); + pair b=(Xmax,y); + pair a2=(Xmin,y2); + pair b2=(Xmax,y2); + + if(finite(a)) { + pic.addPoint(a,min(p)); + pic.addPoint(a,max(p)); + } + + if(finite(b)) { + pic.addPoint(b,min(p)); + pic.addPoint(b,max(p)); + } + + if(finite(a) && finite(b)) { + frame d; + ticks(d,pic.scaling(warn=false),L,side, + (a.x,0)--(b.x,0),(a2.x,0)--(b2.x,0),p,arrow,margin, + ticklocate(a.x,b.x,pic.scale.x),divisor,opposite); + frame f; + if(L.s != "") { + Label L0=L.copy(); + L0.position(0); + add(f,L0); + } + pair pos=a+L.relative()*(b-a); + pic.addBox(pos,pos,(min(f).x,min(d).y),(max(f).x,max(d).y)); + } + } + + // Process any queued y axis bound calculation requests. + for(int i=0; i < pic.scale.y.bound.length; ++i) + pic.scale.y.bound[i](); + + pic.scale.y.bound.delete(); + + bounds(); + + // Request another x bounds calculation before final picture scaling. + pic.scale.x.bound.push(bounds); +} + +// An internal routine to draw a y axis at a particular x value. +void yaxisAt(picture pic=currentpicture, Label L="", axis axis, + real ymin=-infinity, real ymax=infinity, pen p=currentpen, + ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin, + bool above=true, bool opposite=false) +{ + real x=axis.value; + real x2; + Label L=L.copy(); + int[] divisor=copy(axis.ydivisor); + pair side=axis.side; + int type=axis.type; + + pic.add(new void (frame f, transform t, transform T, pair lb, pair rt) { + transform tinv=inverse(t); + pair a=ymin == -infinity ? tinv*(xtrans(t,x),lb.y-min(p).y) : (x,ymin); + pair b=ymax == infinity ? tinv*(xtrans(t,x),rt.y-max(p).y) : (x,ymax); + pair a2=ymin == -infinity ? tinv*(xtrans(t,x2),lb.y-min(p).y) : (x2,ymin); + pair b2=ymax == infinity ? tinv*(xtrans(t,x2),rt.y-max(p).y) : (x2,ymax); + + if(ymin == -infinity || ymax == infinity) { + bounds my=autoscale(a.y,b.y,pic.scale.y.scale); + pic.scale.y.tickMin=my.min; + pic.scale.y.tickMax=my.max; + divisor=my.divisor; + } + + real fuzz=epsilon*max(abs(a.y),abs(b.y)); + a -= (0,fuzz); + b += (0,fuzz); + + frame d; + ticks(d,t,L,side,a--b,finite(x2) ? a2--b2 : nullpath,p,arrow,margin, + ticklocate(a.y,b.y,pic.scale.y),divisor,opposite); + (above ? add : prepend)(f,t*T*tinv*d); + }); + + void bounds() { + if(type == Min) + x=pic.scale.x.automin() ? tickMin(pic).x : pic.userMin.x; + else if(type == Max) + x=pic.scale.x.automax() ? tickMax(pic).x : pic.userMax.x; + else if(type == Both) { + x2=pic.scale.x.automax() ? tickMax(pic).x : pic.userMax.x; + x=opposite ? x2 : + (pic.scale.x.automin() ? tickMin(pic).x : pic.userMin.x); + } + + real Ymin=finite(ymin) ? ymin : pic.userMin.y; + real Ymax=finite(ymax) ? ymax : pic.userMax.y; + + pair a=(x,Ymin); + pair b=(x,Ymax); + pair a2=(x2,Ymin); + pair b2=(x2,Ymax); + + if(finite(a)) { + pic.addPoint(a,min(p)); + pic.addPoint(a,max(p)); + } + + if(finite(b)) { + pic.addPoint(b,min(p)); + pic.addPoint(b,max(p)); + } + + if(finite(a) && finite(b)) { + frame d; + ticks(d,pic.scaling(warn=false),L,side, + (0,a.y)--(0,b.y),(0,a2.y)--(0,b2.y),p,arrow,margin, + ticklocate(a.y,b.y,pic.scale.y),divisor,opposite); + frame f; + if(L.s != "") { + Label L0=L.copy(); + L0.position(0); + add(f,L0); + } + pair pos=a+L.relative()*(b-a); + pic.addBox(pos,pos,(min(d).x,min(f).y),(max(d).x,max(f).y)); + } + } + + // Process any queued x axis bound calculation requests. + for(int i=0; i < pic.scale.x.bound.length; ++i) + pic.scale.x.bound[i](); + + pic.scale.x.bound.delete(); + + bounds(); + + // Request another y bounds calculation before final picture scaling. + pic.scale.y.bound.push(bounds); +} + +// Set the x limits of a picture. +void xlimits(picture pic=currentpicture, real min=-infinity, real max=infinity, + bool crop=NoCrop) +{ + if(min > max) return; + + pic.scale.x.automin=min <= -infinity; + pic.scale.x.automax=max >= infinity; + + bounds mx; + if(pic.scale.x.automin() || pic.scale.x.automax()) + mx=autoscale(pic.userMin.x,pic.userMax.x,pic.scale.x.scale); + + if(pic.scale.x.automin) { + if(pic.scale.x.automin()) pic.userMinx(mx.min); + } else pic.userMinx(pic.scale.x.T(min)); + + if(pic.scale.x.automax) { + if(pic.scale.x.automax()) pic.userMaxx(mx.max); + } else pic.userMaxx(pic.scale.x.T(max)); + + if(crop) { + pair userMin=pic.userMin(); + pair userMax=pic.userMax(); + pic.bounds.xclip(userMin.x,userMax.x); + pic.clip(new void (frame f, transform t, transform T, pair, pair) { + frame Tinvf=T == identity() ? f : t*inverse(T)*inverse(t)*f; + clip(f,T*box(((t*userMin).x,(min(Tinvf)).y), + ((t*userMax).x,(max(Tinvf)).y))); + }); + } +} + +// Set the y limits of a picture. +void ylimits(picture pic=currentpicture, real min=-infinity, real max=infinity, + bool crop=NoCrop) +{ + if(min > max) return; + + pic.scale.y.automin=min <= -infinity; + pic.scale.y.automax=max >= infinity; + + bounds my; + if(pic.scale.y.automin() || pic.scale.y.automax()) + my=autoscale(pic.userMin.y,pic.userMax.y,pic.scale.y.scale); + + if(pic.scale.y.automin) { + if(pic.scale.y.automin()) pic.userMiny(my.min); + } else pic.userMiny(pic.scale.y.T(min)); + + if(pic.scale.y.automax) { + if(pic.scale.y.automax()) pic.userMaxy(my.max); + } else pic.userMaxy(pic.scale.y.T(max)); + + if(crop) { + pair userMin=pic.userMin(); + pair userMax=pic.userMax(); + pic.bounds.yclip(userMin.y,userMax.y); + pic.clip(new void (frame f, transform t, transform T, pair, pair) { + frame Tinvf=T == identity() ? f : t*inverse(T)*inverse(t)*f; + clip(f,T*box(((min(Tinvf)).x,(t*userMin).y), + ((max(Tinvf)).x,(t*userMax).y))); + }); + } +} + +// Crop a picture to the current user-space picture limits. +void crop(picture pic=currentpicture) +{ + xlimits(pic,false); + ylimits(pic,false); + if(pic.userSetx && pic.userSety) + clip(pic,box(pic.userMin(),pic.userMax())); +} + +// Restrict the x and y limits to box(min,max). +void limits(picture pic=currentpicture, pair min, pair max, bool crop=NoCrop) +{ + xlimits(pic,min.x,max.x); + ylimits(pic,min.y,max.y); + if(crop && pic.userSetx && pic.userSety) + clip(pic,box(pic.userMin(),pic.userMax())); +} + +// Internal routine to autoscale the user limits of a picture. +void autoscale(picture pic=currentpicture, axis axis) +{ + if(!pic.scale.set) { + bounds mx,my; + pic.scale.set=true; + + if(pic.userSetx) { + mx=autoscale(pic.userMin.x,pic.userMax.x,pic.scale.x.scale); + if(pic.scale.x.scale.logarithmic && + floor(pic.userMin.x) == floor(pic.userMax.x)) { + if(pic.scale.x.automin()) + pic.userMinx(floor(pic.userMin.x)); + if(pic.scale.x.automax()) + pic.userMaxx(ceil(pic.userMax.x)); + } + } else {mx.min=mx.max=0; pic.scale.set=false;} + + if(pic.userSety) { + my=autoscale(pic.userMin.y,pic.userMax.y,pic.scale.y.scale); + if(pic.scale.y.scale.logarithmic && + floor(pic.userMin.y) == floor(pic.userMax.y)) { + if(pic.scale.y.automin()) + pic.userMiny(floor(pic.userMin.y)); + if(pic.scale.y.automax()) + pic.userMaxy(ceil(pic.userMax.y)); + } + } else {my.min=my.max=0; pic.scale.set=false;} + + pic.scale.x.tickMin=mx.min; + pic.scale.x.tickMax=mx.max; + pic.scale.y.tickMin=my.min; + pic.scale.y.tickMax=my.max; + axis.xdivisor=mx.divisor; + axis.ydivisor=my.divisor; + } +} + +// Draw an x axis. +void xaxis(picture pic=currentpicture, Label L="", axis axis=YZero, + real xmin=-infinity, real xmax=infinity, pen p=currentpen, + ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin, + bool above=false) +{ + if(xmin > xmax) return; + + if(pic.scale.x.automin && xmin > -infinity) pic.scale.x.automin=false; + if(pic.scale.x.automax && xmax < infinity) pic.scale.x.automax=false; + + if(!pic.scale.set) { + axis(pic,axis); + autoscale(pic,axis); + } + + Label L=L.copy(); + bool newticks=false; + + if(xmin != -infinity) { + xmin=pic.scale.x.T(xmin); + newticks=true; + } + + if(xmax != infinity) { + xmax=pic.scale.x.T(xmax); + newticks=true; + } + + if(newticks && pic.userSetx && ticks != NoTicks) { + if(xmin == -infinity) xmin=pic.userMin.x; + if(xmax == infinity) xmax=pic.userMax.x; + bounds mx=autoscale(xmin,xmax,pic.scale.x.scale); + pic.scale.x.tickMin=mx.min; + pic.scale.x.tickMax=mx.max; + axis.xdivisor=mx.divisor; + } + + axis(pic,axis); + + if(xmin == -infinity && !axis.extend) { + if(pic.scale.set) + xmin=pic.scale.x.automin() ? pic.scale.x.tickMin : + max(pic.scale.x.tickMin,pic.userMin.x); + else xmin=pic.userMin.x; + } + + if(xmax == infinity && !axis.extend) { + if(pic.scale.set) + xmax=pic.scale.x.automax() ? pic.scale.x.tickMax : + min(pic.scale.x.tickMax,pic.userMax.x); + else xmax=pic.userMax.x; + } + + if(L.defaultposition) L.position(axis.position); + L.align(L.align,axis.align); + + xaxisAt(pic,L,axis,xmin,xmax,p,ticks,arrow,margin,above); + if(axis.type == Both) + xaxisAt(pic,L,axis,xmin,xmax,p,ticks,arrow,margin,above,true); +} + +// Draw a y axis. +void yaxis(picture pic=currentpicture, Label L="", axis axis=XZero, + real ymin=-infinity, real ymax=infinity, pen p=currentpen, + ticks ticks=NoTicks, arrowbar arrow=None, margin margin=NoMargin, + bool above=false) +{ + if(ymin > ymax) return; + + if(pic.scale.y.automin && ymin > -infinity) pic.scale.y.automin=false; + if(pic.scale.y.automax && ymax < infinity) pic.scale.y.automax=false; + + if(!pic.scale.set) { + axis(pic,axis); + autoscale(pic,axis); + } + + Label L=L.copy(); + bool newticks=false; + + if(ymin != -infinity) { + ymin=pic.scale.y.T(ymin); + newticks=true; + } + + if(ymax != infinity) { + ymax=pic.scale.y.T(ymax); + newticks=true; + } + + if(newticks && pic.userSety && ticks != NoTicks) { + if(ymin == -infinity) ymin=pic.userMin.y; + if(ymax == infinity) ymax=pic.userMax.y; + bounds my=autoscale(ymin,ymax,pic.scale.y.scale); + pic.scale.y.tickMin=my.min; + pic.scale.y.tickMax=my.max; + axis.ydivisor=my.divisor; + } + + axis(pic,axis); + + if(ymin == -infinity && !axis.extend) { + if(pic.scale.set) + ymin=pic.scale.y.automin() ? pic.scale.y.tickMin : + max(pic.scale.y.tickMin,pic.userMin.y); + else ymin=pic.userMin.y; + } + + + if(ymax == infinity && !axis.extend) { + if(pic.scale.set) + ymax=pic.scale.y.automax() ? pic.scale.y.tickMax : + min(pic.scale.y.tickMax,pic.userMax.y); + else ymax=pic.userMax.y; + } + + if(L.defaultposition) L.position(axis.position); + L.align(L.align,axis.align); + + if(L.defaulttransform) { + frame f; + add(f,Label(L.s,(0,0),L.p)); + if(length(max(f)-min(f)) > ylabelwidth*fontsize(L.p)) + L.transform(rotate(90)); + } + + yaxisAt(pic,L,axis,ymin,ymax,p,ticks,arrow,margin,above); + if(axis.type == Both) + yaxisAt(pic,L,axis,ymin,ymax,p,ticks,arrow,margin,above,true); +} + +// Draw x and y axes. +void axes(picture pic=currentpicture, Label xlabel="", Label ylabel="", + pair min=(-infinity,-infinity), pair max=(infinity,infinity), + pen p=currentpen, arrowbar arrow=None, margin margin=NoMargin, + bool above=false) +{ + xaxis(pic,xlabel,min.x,max.x,p,arrow,margin,above); + yaxis(pic,ylabel,min.y,max.y,p,arrow,margin,above); +} + +// Draw a yaxis at x. +void xequals(picture pic=currentpicture, Label L="", real x, + bool extend=false, real ymin=-infinity, real ymax=infinity, + pen p=currentpen, ticks ticks=NoTicks, + arrowbar arrow=None, margin margin=NoMargin, bool above=true) +{ + yaxis(pic,L,XEquals(x,extend),ymin,ymax,p,ticks,arrow,margin,above); +} + +// Draw an xaxis at y. +void yequals(picture pic=currentpicture, Label L="", real y, + bool extend=false, real xmin=-infinity, real xmax=infinity, + pen p=currentpen, ticks ticks=NoTicks, + arrowbar arrow=None, margin margin=NoMargin, bool above=true) +{ + xaxis(pic,L,YEquals(y,extend),xmin,xmax,p,ticks,arrow,margin,above); +} + +pair Scale(picture pic=currentpicture, pair z) +{ + return (pic.scale.x.T(z.x),pic.scale.y.T(z.y)); +} + +real ScaleX(picture pic=currentpicture, real x) +{ + return pic.scale.x.T(x); +} + +real ScaleY(picture pic=currentpicture, real y) +{ + return pic.scale.y.T(y); +} + +// Draw a tick of length size at pair z in direction dir using pen p. +void tick(picture pic=currentpicture, pair z, pair dir, real size=Ticksize, + pen p=currentpen) +{ + pic.add(new void (frame f, transform t) { + pair tz=t*z; + draw(f,tz--tz+unit(dir)*size,p); + }); + pic.addPoint(z,p); + pic.addPoint(z,unit(dir)*size,p); +} + +void xtick(picture pic=currentpicture, explicit pair z, pair dir=N, + real size=Ticksize, pen p=currentpen) +{ + tick(pic,Scale(pic,z),dir,size,p); +} + +void xtick(picture pic=currentpicture, real x, pair dir=N, + real size=Ticksize, pen p=currentpen) +{ + xtick(pic,(x,pic.scale.y.scale.logarithmic ? 1 : 0),dir,size,p); +} + +void ytick(picture pic=currentpicture, explicit pair z, pair dir=E, + real size=Ticksize, pen p=currentpen) +{ + xtick(pic,z,dir,size,p); +} + +void ytick(picture pic=currentpicture, real y, pair dir=E, + real size=Ticksize, pen p=currentpen) +{ + xtick(pic,(pic.scale.x.scale.logarithmic ? 1 : 0,y),dir,size,p); +} + +void tick(picture pic=currentpicture, Label L, real value, explicit pair z, + pair dir, string format="", real size=Ticksize, pen p=currentpen) +{ + Label L=L.copy(); + L.position(z); + L.align(L.align,-dir); + if(shift(L.T)*0 == 0) + L.T=shift(dot(dir,L.align.dir) > 0 ? dir*size : + ticklabelshift(L.align.dir,p))*L.T; + L.p(p); + if(L.s == "") L.s=format(format == "" ? defaultformat : format,value); + L.s=baseline(L.s,baselinetemplate); + add(pic,L); + xtick(pic,z,dir,size,p); +} + +void xtick(picture pic=currentpicture, Label L, explicit pair z, pair dir=N, + string format="", real size=Ticksize, pen p=currentpen) +{ + tick(pic,L,z.x,z,dir,format,size,p); +} + +void xtick(picture pic=currentpicture, Label L, real x, pair dir=N, + string format="", real size=Ticksize, pen p=currentpen) +{ + xtick(pic,L,(x,pic.scale.y.scale.logarithmic ? 1 : 0),dir,size,p); +} + +void ytick(picture pic=currentpicture, Label L, explicit pair z, pair dir=E, + string format="", real size=Ticksize, pen p=currentpen) +{ + tick(pic,L,z.y,z,dir,format,size,p); +} + +void ytick(picture pic=currentpicture, Label L, real y, pair dir=E, + string format="", real size=Ticksize, pen p=currentpen) +{ + ytick(pic,L,(pic.scale.x.scale.logarithmic ? 1 : 0,y),dir,format,size,p); +} + +private void label(picture pic, Label L, pair z, real x, align align, + string format, pen p) +{ + Label L=L.copy(); + L.position(z); + L.align(align); + L.p(p); + if(shift(L.T)*0 == 0) + L.T=shift(ticklabelshift(L.align.dir,L.p))*L.T; + if(L.s == "") L.s=format(format == "" ? defaultformat : format,x); + L.s=baseline(L.s,baselinetemplate); + add(pic,L); +} + +// Put a label on the x axis. +void labelx(picture pic=currentpicture, Label L="", explicit pair z, + align align=S, string format="", pen p=nullpen) +{ + label(pic,L,Scale(pic,z),z.x,align,format,p); +} + +void labelx(picture pic=currentpicture, Label L="", real x, + align align=S, string format="", pen p=nullpen) +{ + labelx(pic,L,(x,pic.scale.y.scale.logarithmic ? 1 : 0),align,format,p); +} + +void labelx(picture pic=currentpicture, Label L, + string format="", explicit pen p=currentpen) +{ + labelx(pic,L,L.position.position,format,p); +} + +// Put a label on the y axis. +void labely(picture pic=currentpicture, Label L="", explicit pair z, + align align=W, string format="", pen p=nullpen) +{ + label(pic,L,Scale(pic,z),z.y,align,format,p); +} + +void labely(picture pic=currentpicture, Label L="", real y, + align align=W, string format="", pen p=nullpen) +{ + labely(pic,L,(pic.scale.x.scale.logarithmic ? 1 : 0,y),align,format,p); +} + +void labely(picture pic=currentpicture, Label L, + string format="", explicit pen p=nullpen) +{ + labely(pic,L,L.position.position,format,p); +} + +private string noprimary="Primary axis must be drawn before secondary axis"; + +// Construct a secondary X axis +picture secondaryX(picture primary=currentpicture, void f(picture)) +{ + if(!primary.scale.set) abort(noprimary); + picture pic; + if(primary.userMax.x == primary.userMin.x) return pic; + f(pic); + if(!pic.userSetx) return pic; + bounds a=autoscale(pic.userMin.x,pic.userMax.x,pic.scale.x.scale); + real bmin=pic.scale.x.automin() ? a.min : pic.userMin.x; + real bmax=pic.scale.x.automax() ? a.max : pic.userMax.x; + + real denom=bmax-bmin; + if(denom != 0) { + pic.erase(); + real m=(primary.userMax.x-primary.userMin.x)/denom; + pic.scale.x.postscale=Linear(m,bmin-primary.userMin.x/m); + pic.scale.set=true; + pic.scale.x.tickMin=pic.scale.x.postscale.T(a.min); + pic.scale.x.tickMax=pic.scale.x.postscale.T(a.max); + pic.scale.y.tickMin=primary.userMin.y; + pic.scale.y.tickMax=primary.userMax.y; + axis.xdivisor=a.divisor; + f(pic); + } + pic.userCopy(primary); + return pic; +} + +// Construct a secondary Y axis +picture secondaryY(picture primary=currentpicture, void f(picture)) +{ + if(!primary.scale.set) abort(noprimary); + picture pic; + if(primary.userMax.y == primary.userMin.y) return pic; + f(pic); + if(!pic.userSety) return pic; + bounds a=autoscale(pic.userMin.y,pic.userMax.y,pic.scale.y.scale); + real bmin=pic.scale.y.automin() ? a.min : pic.userMin.y; + real bmax=pic.scale.y.automax() ? a.max : pic.userMax.y; + + real denom=bmax-bmin; + if(denom != 0) { + pic.erase(); + real m=(primary.userMax.y-primary.userMin.y)/denom; + pic.scale.y.postscale=Linear(m,bmin-primary.userMin.y/m); + pic.scale.set=true; + pic.scale.x.tickMin=primary.userMin.x; + pic.scale.x.tickMax=primary.userMax.x; + pic.scale.y.tickMin=pic.scale.y.postscale.T(a.min); + pic.scale.y.tickMax=pic.scale.y.postscale.T(a.max); + axis.ydivisor=a.divisor; + f(pic); + } + pic.userCopy(primary); + return pic; +} + +typedef guide graph(pair f(real), real, real, int); +typedef guide[] multigraph(pair f(real), real, real, int); + +graph graph(interpolate join) +{ + return new guide(pair f(real), real a, real b, int n) { + real width=b-a; + return n == 0 ? join(f(a)) : + join(...sequence(new guide(int i) {return f(a+(i/n)*width);},n+1)); + }; +} + +multigraph graph(interpolate join, bool3 cond(real)) +{ + return new guide[](pair f(real), real a, real b, int n) { + real width=b-a; + if(n == 0) return new guide[] {join(cond(a) ? f(a) : nullpath)}; + guide[] G; + guide[] g; + for(int i=0; i < n+1; ++i) { + real t=a+(i/n)*width; + bool3 b=cond(t); + if(b) + g.push(f(t)); + else { + if(g.length > 0) { + G.push(join(...g)); + g=new guide[] {}; + } + if(b == default) + g.push(f(t)); + } + } + if(g.length > 0) + G.push(join(...g)); + return G; + }; +} + +guide Straight(... guide[])=operator --; +guide Spline(... guide[])=operator ..; + +interpolate Hermite(splinetype splinetype) +{ + return new guide(... guide[] a) { + int n=a.length; + if(n == 0) return nullpath; + real[] x,y; + guide G; + for(int i=0; i < n; ++i) { + guide g=a[i]; + int m=size(g); + if(m == 0) continue; + pair z=point(g,0); + x.push(z.x); + y.push(z.y); + if(m > 1) { + G=G..hermite(x,y,splinetype) & g; + pair z=point(g,m); + x=new real[] {z.x}; + y=new real[] {z.y}; + } + } + return G & hermite(x,y,splinetype); + }; +} + +interpolate Hermite=Hermite(defaultspline); + +guide graph(picture pic=currentpicture, real f(real), real a, real b, + int n=ngraph, real T(real)=identity, interpolate join=operator --) +{ + if(T == identity) + return graph(join)(new pair(real x) { + return (x,pic.scale.y.T(f(pic.scale.x.Tinv(x))));}, + pic.scale.x.T(a),pic.scale.x.T(b),n); + else + return graph(join)(new pair(real x) { + return Scale(pic,(T(x),f(T(x))));}, + a,b,n); +} + +guide[] graph(picture pic=currentpicture, real f(real), real a, real b, + int n=ngraph, real T(real)=identity, + bool3 cond(real), interpolate join=operator --) +{ + if(T == identity) + return graph(join,cond)(new pair(real x) { + return (x,pic.scale.y.T(f(pic.scale.x.Tinv(x))));}, + pic.scale.x.T(a),pic.scale.x.T(b),n); + else + return graph(join,cond)(new pair(real x) { + return Scale(pic,(T(x),f(T(x))));}, + a,b,n); +} + +guide graph(picture pic=currentpicture, real x(real), real y(real), real a, + real b, int n=ngraph, real T(real)=identity, + interpolate join=operator --) +{ + if(T == identity) + return graph(join)(new pair(real t) {return Scale(pic,(x(t),y(t)));},a,b,n); + else + return graph(join)(new pair(real t) { + return Scale(pic,(x(T(t)),y(T(t)))); + },a,b,n); +} + +guide[] graph(picture pic=currentpicture, real x(real), real y(real), real a, + real b, int n=ngraph, real T(real)=identity, bool3 cond(real), + interpolate join=operator --) +{ + if(T == identity) + return graph(join,cond)(new pair(real t) {return Scale(pic,(x(t),y(t)));}, + a,b,n); + else + return graph(join,cond)(new pair(real t) { + return Scale(pic,(x(T(t)),y(T(t))));}, + a,b,n); +} + +guide graph(picture pic=currentpicture, pair z(real), real a, real b, + int n=ngraph, real T(real)=identity, interpolate join=operator --) +{ + if(T == identity) + return graph(join)(new pair(real t) {return Scale(pic,z(t));},a,b,n); + else + return graph(join)(new pair(real t) { + return Scale(pic,z(T(t))); + },a,b,n); +} + +guide[] graph(picture pic=currentpicture, pair z(real), real a, real b, + int n=ngraph, real T(real)=identity, bool3 cond(real), + interpolate join=operator --) +{ + if(T == identity) + return graph(join,cond)(new pair(real t) {return Scale(pic,z(t));},a,b,n); + else + return graph(join,cond)(new pair(real t) { + return Scale(pic,z(T(t))); + },a,b,n); +} + +string differentlengths="attempt to graph arrays of different lengths"; +string conditionlength="condition array has different length than data"; + +void checklengths(int x, int y, string text=differentlengths) +{ + if(x != y) + abort(text+": "+string(x)+" != "+string(y)); +} + +void checkconditionlength(int x, int y) +{ + checklengths(x,y,conditionlength); +} + +guide graph(picture pic=currentpicture, pair[] z, interpolate join=operator --) +{ + int i=0; + return graph(join)(new pair(real) { + pair w=Scale(pic,z[i]); + ++i; + return w; + },0,0,z.length-1); +} + +guide[] graph(picture pic=currentpicture, pair[] z, bool3[] cond, + interpolate join=operator --) +{ + int n=z.length; + int i=0; + pair w; + checkconditionlength(cond.length,n); + bool3 condition(real) { + bool3 b=cond[i]; + if(b != false) w=Scale(pic,z[i]); + ++i; + return b; + } + return graph(join,condition)(new pair(real) {return w;},0,0,n-1); +} + +guide graph(picture pic=currentpicture, real[] x, real[] y, + interpolate join=operator --) +{ + int n=x.length; + checklengths(n,y.length); + int i=0; + return graph(join)(new pair(real) { + pair w=Scale(pic,(x[i],y[i])); + ++i; + return w; + },0,0,n-1); +} + +guide[] graph(picture pic=currentpicture, real[] x, real[] y, bool3[] cond, + interpolate join=operator --) +{ + int n=x.length; + checklengths(n,y.length); + int i=0; + pair w; + checkconditionlength(cond.length,n); + bool3 condition(real) { + bool3 b=cond[i]; + if(b != false) w=Scale(pic,(x[i],y[i])); + ++i; + return b; + } + return graph(join,condition)(new pair(real) {return w;},0,0,n-1); +} + +// Connect points in z into segments corresponding to consecutive true elements +// of b using interpolation operator join. +path[] segment(pair[] z, bool[] cond, interpolate join=operator --) +{ + checkconditionlength(cond.length,z.length); + int[][] segment=segment(cond); + return sequence(new path(int i) {return join(...z[segment[i]]);}, + segment.length); +} + +pair polar(real r, real theta) +{ + return r*expi(theta); +} + +guide polargraph(picture pic=currentpicture, real r(real), real a, real b, + int n=ngraph, interpolate join=operator --) +{ + return graph(join)(new pair(real theta) { + return Scale(pic,polar(r(theta),theta)); + },a,b,n); +} + +void errorbar(picture pic, pair z, pair dp, pair dm, pen p=currentpen, + real size=0) +{ + real dmx=-abs(dm.x); + real dmy=-abs(dm.y); + real dpx=abs(dp.x); + real dpy=abs(dp.y); + if(dmx != dpx) draw(pic,Scale(pic,z+(dmx,0))--Scale(pic,z+(dpx,0)),p, + Bars(size)); + if(dmy != dpy) draw(pic,Scale(pic,z+(0,dmy))--Scale(pic,z+(0,dpy)),p, + Bars(size)); +} + +void errorbars(picture pic=currentpicture, pair[] z, pair[] dp, pair[] dm={}, + bool[] cond={}, pen p=currentpen, real size=0) +{ + if(dm.length == 0) dm=dp; + int n=z.length; + checklengths(n,dm.length); + checklengths(n,dp.length); + bool all=cond.length == 0; + if(!all) + checkconditionlength(cond.length,n); + for(int i=0; i < n; ++i) { + if(all || cond[i]) + errorbar(pic,z[i],dp[i],dm[i],p,size); + } +} + +void errorbars(picture pic=currentpicture, real[] x, real[] y, + real[] dpx, real[] dpy, real[] dmx={}, real[] dmy={}, + bool[] cond={}, pen p=currentpen, real size=0) +{ + if(dmx.length == 0) dmx=dpx; + if(dmy.length == 0) dmy=dpy; + int n=x.length; + checklengths(n,y.length); + checklengths(n,dpx.length); + checklengths(n,dpy.length); + checklengths(n,dmx.length); + checklengths(n,dmy.length); + bool all=cond.length == 0; + if(!all) + checkconditionlength(cond.length,n); + for(int i=0; i < n; ++i) { + if(all || cond[i]) + errorbar(pic,(x[i],y[i]),(dpx[i],dpy[i]),(dmx[i],dmy[i]),p,size); + } +} + +void errorbars(picture pic=currentpicture, real[] x, real[] y, + real[] dpy, bool[] cond={}, pen p=currentpen, real size=0) +{ + errorbars(pic,x,y,0*x,dpy,cond,p,size); +} + +// Return a vector field on path g, specifying the vector as a function of the +// relative position along path g in [0,1]. +picture vectorfield(path vector(real), path g, int n, bool truesize=false, + pen p=currentpen, arrowbar arrow=Arrow, + margin margin=PenMargin) +{ + picture pic; + for(int i=0; i < n; ++i) { + real x=(n == 1) ? 0.5 : i/(n-1); + if(truesize) + draw(relpoint(g,x),pic,vector(x),p,arrow); + else + draw(pic,shift(relpoint(g,x))*vector(x),p,arrow,margin); + } + return pic; +} + +real maxlength(pair a, pair b, int nx, int ny) +{ + return min((b.x-a.x)/nx,(b.y-a.y)/ny); +} + +// return a vector field over box(a,b). +picture vectorfield(path vector(pair), pair a, pair b, + int nx=nmesh, int ny=nx, bool truesize=false, + real maxlength=truesize ? 0 : maxlength(a,b,nx,ny), + bool cond(pair z)=null, pen p=currentpen, + arrowbar arrow=Arrow, margin margin=PenMargin) +{ + picture pic; + real dx=1/nx; + real dy=1/ny; + bool all=cond == null; + real scale; + + if(maxlength > 0) { + real size(pair z) { + path g=vector(z); + return abs(point(g,size(g)-1)-point(g,0)); + } + real max=size((0,0)); + for(int i=0; i <= nx; ++i) { + real x=interp(a.x,b.x,i*dx); + for(int j=0; j <= ny; ++j) + max=max(max,size((x,interp(a.y,b.y,j*dy)))); + } + scale=max > 0 ? maxlength/max : 1; + } else scale=1; + + for(int i=0; i <= nx; ++i) { + real x=interp(a.x,b.x,i*dx); + for(int j=0; j <= ny; ++j) { + real y=interp(a.y,b.y,j*dy); + pair z=(x,y); + if(all || cond(z)) { + path g=scale(scale)*vector(z); + if(truesize) + draw(z,pic,g,p,arrow); + else + draw(pic,shift(z)*g,p,arrow,margin); + } + } + } + return pic; +} + +// True arc +path Arc(pair c, real r, real angle1, real angle2, bool direction, + int n=nCircle) +{ + angle1=radians(angle1); + angle2=radians(angle2); + if(angle1 >= angle2 && direction) angle1 -= 2pi; + if(angle2 >= angle1 && !direction) angle2 -= 2pi; + return shift(c)*polargraph(new real(real t){return r;},angle1,angle2,n, + operator ..); +} + +path Arc(pair c, real r, real angle1, real angle2, int n=nCircle) +{ + return Arc(c,r,angle1,angle2,angle2 >= angle1 ? CCW : CW,n); +} + +path Arc(pair c, explicit pair z1, explicit pair z2, bool direction=CCW, + int n=nCircle) +{ + return Arc(c,abs(z1-c),degrees(z1-c),degrees(z2-c),direction,n); +} + +// True circle +path Circle(pair c, real r, int n=nCircle) +{ + return Arc(c,r,0,360,n)&cycle; +} |