% sriyantraquest.mp % L. Nobre G. % 2009 prologues := 3; def ap( expr poia, poib ) = begingroup pair dirv; dirv = 2*u*unitvector(poib-poia); ( (poia-dirv)--(poia+dirv) ) endgroup enddef; def masstri( expr A, B, C ) = begingroup numeric mab, mbc, mca, mto; mab = 0.005*abs( A - B ); mbc = 0.005*abs( B - C ); mca = 0.005*abs( C - A ); mto = mab+mbc+mca; ( ( 0.5*ypart (mab[A,B]+mbc[B,C]+mca[C,A]), mto ) ) endgroup enddef; beginfig(1); numeric u, gp[], i, j, k, p[], g[], costval, mincost, marg, massto; numeric deviatval, reducval, concurrpoiscale, auxsum, counter, lastmin; numeric mindist, deviatmax; path outercircle; pair tp, tq; boolean ba, bb, bc, bd, doneflag; outercircle = fullcircle scaled 2; u = 5cm; mindist = 0.1; deviatmax = 0.015; concurrpoiscale = 3pt; counter = 1; doneflag = false; % deviatval = 0.5; % reducval = 0.75; % marg = 0.2; % costval = 0.0043; % mincost = 0.00429; % gp1 = 0.76746; % gp2 = 0.34315; % gp3 = -0.22505; % gp4 = 0.24968; % gp5 = 0.5209; % deviatval = 0.525; % reducval = 0.75; % marg = 0.2; % costval = 0.00273; % mincost = 0.00272; % gp1 = 0.78983; % gp2 = 0.4125; % gp3 = -0.27725; % gp4 = 0.16902; % gp5 = 0.44658; % deviatval = 0.55; % reducval = 0.75; % marg = 0.2; % costval = 0.0048; % mincost = 0.0047; % gp1 = 0.70071; % gp2 = 0.3247; % gp3 = -0.20796; % gp4 = 0.24506; % gp5 = 0.5587; % deviatval = 0.575; % reducval = 0.75; % marg = 0.05; % costval = 0.0004; % mincost = 0.0003; % gp1 = 0.82916; % gp2 = 0.41296; % gp3 = -0.25703; % gp4 = 0.19568; % gp5 = 0.44472; % deviatval = 0.6; % reducval = 0.5; % marg = 0.025; % costval = 0.00636; % mincost = 0.00635; % gp1 = 0.76761; % gp2 = 0.2831; % gp3 = -0.15158; % gp4 = 0.42883; % gp5 = 0.60568; % deviatval = 0.625; % reducval = 0.5; % marg = 0.05; % costval = 0.0004; % mincost = 0.0003; % gp1 = 0.83; % gp2 = 0.38884; % gp3 = -0.23164; % gp4 = 0.23488; % gp5 = 0.47797; % deviatval = 0.65; % reducval = 0.9; % marg = 0.01; % costval = 0.009; % mincost = 0.008; % gp1 = 0.74031; % gp2 = 0.2599; % gp3 = -0.12587; % gp4 = 0.50862; % gp5 = 0.64258; % deviatval = 0.7; % reducval = 0.75; % marg = 0.2; % costval = 0.00292; % mincost = 0.00291; % gp1 = 0.81487; % gp2 = 0.34457; % gp3 = -0.17764; % gp4 = 0.32793; % gp5 = 0.53682; % deviatval = 0.75; % reducval = 0.75; % marg = 0.2; % costval = 0.01; % mincost = 0.0097; % gp1 = 0.71198; % gp2 = 0.23265; % gp3 = -0.10187; % gp4 = 0.61037; % gp5 = 0.69113; % deviatval = 0.8; % reducval = 0.9; % marg = 0.25; % costval = 0.00038; % mincost = 0.00037; % gp1 = 0.79955; % gp2 = 0.28337; % gp3 = -0.12138; % gp4 = 0.53542; % gp5 = 0.6276; deviatval = 0.8; reducval = 0.75; marg = 0.005; costval = 0.0079; mincost = 0.0078; gp1 = 0.65149; gp2 = 0.20883; gp3 = -0.09235; gp4 = 0.65298; gp5 = 0.74846; % deviatval = 0.85; % reducval = 0.75; % marg = 0.015; % costval = 0.0005; % mincost = 0.00049; % gp1 = 0.76701; % gp2 = 0.28557; % gp3 = -0.12613; % gp4 = 0.48254; % gp5 = 0.63602; % deviatval = 0.85; % reducval = 0.75; % marg = 0.025; % costval = 0.00213; % mincost = 0.00212; % gp1 = 0.6715; % gp2 = 0.21352; % gp3 = -0.08372; % gp4 = 0.69928; % gp5 = 0.73703; % deviatval = 0.9; % reducval = 0.5; % marg = 0.005; % costval = 0.0019; % mincost = 0.00189; % gp1 = 0.66177; % gp2 = 0.21542; % gp3 = -0.08813; % gp4 = 0.65677; % gp5 = 0.74582; % deviatval = 0.9; % reducval = 0.5; % marg = 0.05; % costval = 0.00033; % mincost = 0.00032; % gp1 = 0.80038; % gp2 = 0.37895; % gp3 = -0.18217; % gp4 = 0.26828; % gp5 = 0.51485; lastmin = mincost+0.0001; for j=1 upto 5: g[j] := gp[j]; endfor; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=2000 step -400 until 400: for k=1 upto i: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % forever: % for k=1 upto 256: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% p1 := gp1+uniformdeviate(marg)-0.5*marg; p2 := gp2+uniformdeviate(marg)-0.5*marg; p4 := gp4+uniformdeviate(marg)-0.5*marg; p5 := gp5+uniformdeviate(marg)-0.5*marg; p3 := gp3+uniformdeviate(marg)-0.5*marg; % p3 := p1-sqrt(3)*p4+uniformdeviate(marg)-0.5*marg; % p3 := p1-sqrt(3)*p4; if p3 < 0: pair q[]; q1 = (p4,p1); q2 = (-p4,p1); q3 = (0,p3); q10 = ((-1,-p2)--(1,-p2)) intersectionpoint reverse outercircle; q11 = (-xpart q10,-p2); q12 = (0,1); path a, b; a = ap(q1,q3); b = ap(q10,q12); tp := a intersectiontimes b; if tp <> (-1,-1): q13 = a intersectionpoint b; q17 = (-xpart q13,ypart q13); q24 = 0.5[q13,q17]; q30 = 0.5[q10,q11]; q16 = (0,-p1); ypart q15 = ypart q13; xpart q15 = p5; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% q14 = whatever[q15,q16]; q14 = whatever[q10,q12]; path a; a = (-1,ypart q14)--(1,ypart q14); tp := a intersectiontimes outercircle; if tp <> (-1,-1): q4 = a intersectionpoint outercircle; q5 = (-xpart q4,ypart q14); q6 = (0,-1); q18 = (-xpart q15,ypart q15); q36 = whatever[q1,q3]; q36 = whatever[q4,q5]; q35 = whatever[q16,q15]; q35 = whatever[q10,q11]; q9 = 0.5[q4,q5]; q7 = whatever[q9,q35]; ypart q7 = ypart q16; q8 = (-xpart q7,ypart q7); q37 = whatever[q4,q6]; q37 = whatever[q10,q11]; q38 = whatever[q7,q9]; q38 = whatever[q4,q6]; xpart q19 = 0; q19 = whatever[q37,q36]; ypart q40 = ypart q19; q40 = whatever[q3,q13]; q39 = (-xpart q40,ypart q40); ypart q20 = ypart q38; q20 = whatever[q37,q19]; q21 = (-xpart q20,ypart q20); path a, b; a = q3--(1,p3); b = ap(q14,q16); tp := a intersectiontimes b; if tp <> (-1,-1): q22 = a intersectionpoint b; q23 = (-xpart q22,ypart q22); q25 = 0.5[q20,q21]; path a, b; a = ap(q3,q1); b = ap(q22,q24); tp := a intersectiontimes b; if tp <> (-1,-1): q26 = a intersectionpoint b; q27 = (-xpart q26,ypart q26); q28 = whatever[q19,q20]; q28 = whatever[q26,q27]; q29 = (-xpart q28,ypart q28); path a, b; a = ap(q7,q9); b = ap(q1,q3); tp := a intersectiontimes b; if tp <> (-1,-1): q31 = a intersectionpoint b; q32 = (-xpart q31,ypart q31); q33 = whatever[q31,q32]; q33 = whatever[q22,q24]; tp := (q25--q28) intersectiontimes (q30--q33); tq := (q25--q28) intersectiontimes (q15--q16); if (tp=(-1,-1)) and (tq=(-1,-1)): path a, b; a = ap(q22,q26); b = ap(q25,q28); tp := a intersectiontimes b; if tp <> (-1,-1): q42 = a intersectionpoint b; if abs(q33-deviatval[q26,q42]) < deviatmax: q34 = (-xpart q33,ypart q33); q41 = whatever*dir(angle(q31-q3)-90); q41 = whatever[q3,q31]; path pe, pf, pg; pe = ap(q7,q9); pf = ap(q28,q25); pg = ap(q23,(1,ypart q22)); pair ta, tb, tc; ta = pe intersectiontimes pf; tb = pf intersectiontimes pg; tc = pg intersectiontimes pe; if (ta<>(-1,-1)) and (tb<>(-1,-1)) and (tc<>(-1,-1)): q100 = pe intersectionpoint pf; q101 = pf intersectionpoint pg; q102 = pg intersectionpoint pe; costval := 0; massto := 0; q999 := masstri(q10,q11,q12); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q6,q5,q4); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q21,q20,q19); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q15,q16,q18); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q7,q8,q9); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q40,q39,q3); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q22,q23,q24); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q29,q28,q25); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(q30,q33,q34); costval := costval+xpart q999; massto := massto + ypart q999; costval := costval/massto; costval := 1+0.08*abs(costval)**2; auxsum := abs(q100-q101)**2+abs(q101-q102)**2; auxsum := 2*(auxsum + abs(q102-q100)**2); costval := costval*(1+auxsum); costval := costval*(1+10*abs(abs(q41)-abs(0.5[q31,q32]))**2); costval := costval*(1+10*abs(q33-deviatval[q26,q42])**2); costval := costval-1; if abs(q40)>(1-mindist): costval := costval+5*(abs(q40)-(1-mindist))**2; fi; if abs(q7)>(1-mindist): costval := costval+5*(abs(q7)-(1-mindist))**2; fi; if abs(q20)>(1-mindist): costval := costval+5*(abs(q20)-(1-mindist))**2; fi; if abs(q25-q16) 75 ); % bc := not ( doneflag or ( lastmin <> mincost ) ); % bd := ( marg < 0.00006 ); % exitif ba or bb or bc or bd; % lastmin := mincost; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% doneflag := false; endfor; show (mincost,angle((gp1-gp3,gp4))); numeric pA, pB, pC, pD, pE; path outercircle; outercircle = fullcircle scaled (2*u); draw outercircle; pA = gp1; pB = gp2; pC = gp3; pD = gp4; pE = gp5; message " gp1 = " & decimal(pA) & "; "; message " gp2 = " & decimal(pB) & "; "; message " gp3 = " & decimal(pC) & "; "; message " gp4 = " & decimal(pD) & "; "; message " gp5 = " & decimal(pE) & "; "; pA := pA*u; pB := pB*u; pC := pC*u; pD := pD*u; pE := pE*u; z1 = (pD,pA); z2 = (-pD,pA); z3 = (0,pC); z10 = ((-u,-pB)--(u,-pB)) intersectionpoint reverse outercircle; z11 = (-x10,-pB); z12 = (0,u); draw z10--z11--z12--cycle; z13 = (z1--z3) intersectionpoint (z10--z12); z17 = (-x13,y13); z24 = 0.5[z13,z17]; z30 = 0.5[z10,z11]; z16 = (0,-pA); y15 = y13; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x15 = pE; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z14 = whatever[z15,z16]; z14 = whatever[z10,z12]; z4 = ((-u,y14)--(u,y14)) intersectionpoint outercircle; z5 = (-x4,y14); z6 = (0,-u); draw z4--z5--z6--cycle; z18 = (-x15,y15); draw z15--z16--z18--cycle; z36 = whatever[z1,z3]; z36 = whatever[z4,z5]; z35 = whatever[z16,z15]; z35 = whatever[z10,z11]; z9 = 0.5[z4,z5]; z7 = whatever[z9,z35]; y7 = y16; z8 = (-x7,y7); draw z7--z8--z9--cycle; z37 = whatever[z4,z6]; z37 = whatever[z10,z11]; z38 = whatever[z7,z9]; z38 = whatever[z4,z6]; x19 = 0; z19 = whatever[z37,z36]; y40 = y19; z40 = whatever[z3,z13]; z39 = (-x40,y40); draw z39--z40--z3--cycle; y20 = y38; z20 = whatever[z37,z19]; z21 = (-x20,y20); draw z19--z20--z21--cycle; z22 = (z3--(u,pC)) intersectionpoint (z14--z16); z23 = (-x22,y22); draw z22--z23--z24--cycle; z25 = 0.5[z20,z21]; z26 = (z3--z1) intersectionpoint (z22--z24); z27 = (-x26,y26); z28 = whatever[z19,z20]; z28 = whatever[z26,z27]; z29 = (-x28,y28); draw z25--z28--z29--cycle; z31 = (z7--z9) intersectionpoint (z1--z3); z32 = (-x31,y31); z33 = whatever[z31,z32]; z33 = whatever[z22,z24]; z34 = (-x33,y33); draw z30--z33--z34--cycle; dotlabels.top(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40); path pe, pf, pg; pe = ap(z7,z9); pf = ap(z28,z25); pg = ap(z23,(u,y22)); filldraw buildcycle( pe, pf, pg ) withcolor red; draw origin withpen pencircle scaled concurrpoiscale withcolor red; draw fullcircle scaled (2*y31) withcolor blue; costval := 0; massto := 0; q999 := masstri(z10,z11,z12); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z6,z5,z4); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z21,z20,z19); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z15,z16,z18); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z7,z8,z9); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z40,z39,z3); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z22,z23,z24); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z29,z28,z25); costval := costval+xpart q999; massto := massto + ypart q999; q999 := masstri(z30,z33,z34); costval := costval+xpart q999; massto := massto + ypart q999; costval := costval/massto; draw (0,costval) withpen pencircle scaled concurrpoiscale withcolor green; endfig; end.