// betatron resonance
// loworder, highorder
// lownux, highnux, lownuy, highnuy
//

// m*nux+n*nuy=k
//

int loworder=1;
int highorder=2;

picture resonancelines(int order, real[] nux, real[] nuy, pen p) {
    picture pic;
    // m*nux+n*nuy = k
    int m,n,i,j;
    int kmin, kmax;
    pair[] c = new pair[4];
    for(int m=-order; m < order; ++m) {
	for (int j = 0; j < 2; ++j) {
	    n = ((-1)^j)*(order-abs(m));
	    write("m,n = ",m,n);
	    //i = abs(m);
	    kmin = -floor(abs(m)*nux[1]+abs(n)*nuy[1]);
	    kmax = floor(abs(m)*nux[1]+abs(n)*nuy[1]);

	    for(int k = kmin; k < kmax; ++k) {
		if(m ==0 && n == 0) break;
		else if (m == 0) {
		    draw((nux[0], k*1.0/n)--(nux[1], k*1.0/n), p);
		    continue;
		}
		else if (n == 0) {
		    draw((k*1.0/m, nuy[0])--(k*1.0/m, nuy[1]), p);
		    continue;
		}		    
		c[0] = (nux[0], (k-nux[0]*m)/n);
		c[1] = (nux[1], (k-nux[1]*m)/n);
		c[2] = ((k-nuy[0]*n)/m, nuy[0]);
		c[3] = ((k-nuy[1]*n)/m, nuy[1]);

		if (c[0].x < nux[0] && c[1].x < nux[0] && c[2].x < nux[0] && c[3].x < nux[0]) continue;
		if (c[0].x > nux[1] && c[1].x > nux[1] && c[2].x > nux[1] && c[3].x > nux[1]) continue;
		if (c[0].y < nuy[0] && c[1].y < nuy[0] && c[2].y < nuy[0] && c[3].y < nuy[0]) continue;
		if (c[0].y > nuy[1] && c[1].y > nuy[1] && c[2].y > nuy[1] && c[3].y > nuy[1]) continue;

		draw(c[0]--c[1]--c[2]--c[3]--cycle, p);
	        //dot(c,4pt+blue);
	    }
	}
    }
    return pic;
}

real[] nux={10.0, 12.0};
//real[] nux={0.0, 2.0};
real[] nuy={0.0, 2.0};

size(600,400);
unitsize(5cm);


draw((nux[0],nuy[0])--(nux[1],nuy[0])--(nux[1],nuy[1])--(nux[0],nuy[1])--cycle,black+2pt);
resonancelines(1, nux, nuy, 1pt+black);
resonancelines(2, nux, nuy, 1pt+black);
resonancelines(3, nux, nuy, 1pt+blue+dotted);
clip((nux[0],nuy[0])--(nux[1],nuy[0])--(nux[1],nuy[1])--(nux[0],nuy[1])--cycle);

for (real i = nux[0]; i <= nux[1]; i=i+1)
    label((string)i, (i,nuy[0]), S);

for (real i = nuy[0]; i <= nuy[1]; i=i+1)
    label((string)i, (nux[0],i), W);