/* built with Studio Sketchpad: 
 *   https://sketchpad.cc
 * 
 * observe the evolution of this sketch: 
 *   https://studio.sketchpad.cc/sp/pad/view/ro.Hbaqzm5iMdE/rev.1839
 * 
 * authors: 
 *   Etienne Balit
 *   

 * license (unless otherwise specified): 
 *   creative commons attribution-share alike 3.0 license.
 *   https://creativecommons.org/licenses/by-sa/3.0/ 
 */ 



void setup() {
    background('white');
    size(600, 600); 
    smooth();
    frameRate(60);
    strokeWeight(1);
    noFill();

    ellipseMode(CENTER);
} 

void draw() {
    faisceau(mouseX, mouseY);
}

void mouseClicked() {
    faisceau(mouseX, mouseY);
}

void faisceau(int ix, int iy) {
    background('white');
    float cibleX = width/2;
    float cibleY = height/2;
    float x = ix;
    float y = iy;
    
    float biss = PI/2;
    float phi = biss;
    float init = biss;
    
    for(i=0; i < 2000; i++){
        stroke(0);
        float c = 20;
        
        float dist1 = dist(cibleX - c, cibleY, x, y);
        float dist2 = dist(cibleX + c, cibleY, x, y);
        
        float dx1 = (cibleX  - c) - x;
        float dx2 = (cibleX + c) - x;
        float dy = cibleY - y;
        
        float bearing1 = atan2(dy, dx1);
        float bearing2 = atan2(dy, dx2);
        
        float lastphi = phi;
        float biss = (bearing1 + bearing2) / 2;
        
        float a = (dist1 + dist2) / 2;
        float b = sqrt(a*a - c*c);
        
        float alpha = dist1 - dist2;
        
        if(i % 200 == 0) {
            stroke(255, 0, 0);
            line( x, y, x - alpha*cos(biss-PI/2), y - alpha*sin(biss-PI/2) );

            stroke(255, 0, 255, 75);
            line(x, y, cibleX  - c, cibleY );
            line(x, y, cibleX + c, cibleY);
            
            stroke(0, 0, 255, 75);
            ellipse(cibleX, cibleY, 2*a, 2*b);
        }
        
        if(i == 0) {
            stroke(0, 0, 0, 75);
            line( x, y, x+1000*cos(biss), y+1000*sin(biss) );  
        }
        
        stroke(0, 255, 0);
        line(x, y, x, y);
        
        float v = 1 - 1/(1 + exp(-0.05*(a-100)));

        float dX = cos(biss) - 0.1*v*alpha*cos(biss-PI/2);
        float dY = sin(biss) - 0.1*v*alpha*sin(biss-PI/2);
        float phi = atan2(dY, dX);
        
        float diff = 10*abs(lastphi-phi);
        float t = (1-cos(diff/2))/2; //1 - 1/(1 + exp(-0.01*(diff-PI)));
        
        x = x + 0.5 * ( (1-t)*cos(phi) + t*cos(lastphi) );
        y = y + 0.5 * ( (1-t)*sin(phi) + t*sin(lastphi) );
        
        float epsilon = dist(x, y, width/2, height/2);
        if(epsilon < 1) { break; }
    }
}