/***********************************************/
/*   Galaxy-Galaxy Collision Simulation        */
/*                                             */
/*             1997.9.26   by Ryoji Matsumoto  */
/***********************************************/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define N 400

main()
{
	double x[1000],y[1000],vx[1000],vy[1000],t,dt,tmax,r,r3,th;
        double k0[4],k1[4],k2[4],k3[4];
        double x1,y1,vx1,vy1,x2,y2,vx2,vy2;
        double x3,y3,vx3,vy3,x4,y4,vx4,vy4;
        double xa,ya,xb,yb,ra,rb,ra3,rb3,mb;
        double vxa,vya,vxb,vyb,xan,yan,vxan,vyan,xbn,ybn,vxbn,vybn;
        int i;

        xa=0.0;
        ya=0.0;
        vxa=0.0;
        vya=0.0;
        xb =2.0;
        yb =-2.0;
        mb =0.3;
	scanf("%lf %lf %lf %lf",&dt,&tmax,&vxb,&vyb);

        for(i=0;i<N/2;i++) {
                r = 0.08*(i/20)+0.2;
                th= 0.31416*(i%20);
                x[i] = xa+r*cos(th);
                y[i] = ya+r*sin(th);
                vx[i]= vxa-sin(th)*sqrt(1.0/r);
                vy[i]= vya+cos(th)*sqrt(1.0/r);
	      }
        for(i=N/2;i<N;i++) {
                r = 0.04*((i-N/2)/20)+0.2;
                th= 0.31416*((i-N/2)%20);
                x[i] = xb+r*cos(th);
                y[i] = yb+r*sin(th);
                vx[i]= vxb-sin(th)*sqrt(mb/r);
                vy[i]= vyb+cos(th)*sqrt(mb/r);
	      }

	for(t=0;t<tmax;t+=dt) {
                rb = sqrt((xa-xb)*(xa-xb)+(ya-yb)*(ya-yb));
                rb3= rb*rb*rb;
                k0[0]= dt*vxa;
                k0[1]= dt*vya;
                k0[2]= dt*(-mb*(xa-xb)/rb3);
                k0[3]= dt*(-mb*(ya-yb)/rb3);
                x1 = xa+k0[0]/2.0;
                y1 = ya+k0[1]/2.0;
                vx1=vxa+k0[2]/2.0;
                vy1=vya+k0[3]/2.0;
                rb = sqrt((x1-xb)*(x1-xb)+(y1-yb)*(y1-yb));
                rb3= rb*rb*rb;
                k1[0]= dt*vx1;
                k1[1]= dt*vy1;
                k1[2]= dt*(-mb*(x1-xb)/rb3);
                k1[3]= dt*(-mb*(y1-yb)/rb3);
                x2 = xa+k1[0]/2.0;
                y2 = ya+k1[1]/2.0;
                vx2=vxa+k1[2]/2.0;
                vy2=vya+k1[3]/2.0;
                rb = sqrt((x2-xb)*(x2-xb)+(y2-yb)*(y2-yb));
                rb3= rb*rb*rb;
                k2[0]= dt*vx2;
                k2[1]= dt*vy2;
                k2[2]= dt*(-mb*(x2-xb)/rb3);
                k2[3]= dt*(-mb*(y2-yb)/rb3);
                x3 = xa+k2[0];
                y3 = ya+k2[1];
                vx3=vxa+k2[2];
                vy3=vya+k2[3];
                rb = sqrt((x3-xb)*(x3-xb)+(y3-yb)*(y3-yb));
                rb3= rb*rb*rb;
                k3[0]= dt*vx3;
                k3[1]= dt*vy3;
                k3[2]= dt*(-mb*(x3-xb)/rb3);
                k3[3]= dt*(-mb*(y3-yb)/rb3);
                xan =xa +(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                yan =ya +(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                vxan=vxa+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                vyan=vya+(k0[3]+2.0*k1[3]+2.0*k2[3]+k3[3])/6.0;

                ra = sqrt((xb-xa)*(xb-xa)+(yb-ya)*(yb-ya));
                ra3= ra*ra*ra;
                k0[0]=dt*vxb;
                k0[1]=dt*vyb;
                k0[2]=dt*(-(xb-xa)/ra3);
                k0[3]=dt*(-(yb-ya)/ra3);
                x1 =xb+k0[0]/2.0;
                y1 =yb+k0[1]/2.0;
                vx1=vxb+k0[2]/2.0;
                vy1=vyb+k0[3]/2.0;
                ra =sqrt((x1-xa)*(x1-xa)+(y1-ya)*(y1-ya));
                ra3= ra*ra*ra;
                k1[0]=dt*vx1;
                k1[1]=dt*vy1;
                k1[2]=dt*(-(x1-xa)/ra3);
                k1[3]=dt*(-(y1-ya)/ra3);
                x2 =xb+k1[0]/2.0;
                y2 =yb+k1[1]/2.0;
                vx2=vxb+k1[2]/2.0;
                vy2=vyb+k1[3]/2.0;
                ra =sqrt((x2-xa)*(x2-xa)+(y2-ya)*(y2-ya));
                ra3=ra*ra*ra;
                k2[0]=dt*vx2;
                k2[1]=dt*vy2;
                k2[2]=dt*(-(x2-xa)/ra3);
                k2[3]=dt*(-(y2-ya)/ra3);
                x3 =xb+k2[0];
                y3 =yb+k2[1];
                vx3=vxb+k2[2];
                vy3=vyb+k2[3];
                ra =sqrt((x3-xa)*(x3-xa)+(y3-ya)*(y3-ya));
                ra3=ra*ra*ra;
                k3[0]=dt*vx3;
                k3[1]=dt*vy3;
                k3[2]=dt*(-(x3-xa)/ra3);
                k3[3]=dt*(-(y3-ya)/ra3);
                xbn =xb +(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                ybn =yb +(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
                vxbn=vxb+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
                vybn=vyb+(k0[3]+2.0*k1[3]+2.0*k2[3]+k3[3])/6.0;

                for(i=0;i<N-1;i++){
                ra = sqrt((x[i]-xa)*(x[i]-xa)+(y[i]-ya)*(y[i]-ya));
                rb = sqrt((x[i]-xb)*(x[i]-xb)+(y[i]-yb)*(y[i]-yb));   
                ra3= ra*ra*ra;
                rb3= rb*rb*rb;
                k0[0]=dt*vx[i];
                k0[1]=dt*vy[i];
                k0[2]=dt*(-(x[i]-xa)/ra3-mb*(x[i]-xb)/rb3);
                k0[3]=dt*(-(y[i]-ya)/ra3-mb*(y[i]-yb)/rb3);
                x1 = x[i]+k0[0]/2.0;
                y1 = y[i]+k0[1]/2.0;
                vx1= vx[i]+k0[2]/2.0;
                vy1= vy[i]+k0[3]/2.0;
                ra = sqrt((x1-xa)*(x1-xa)+(y1-ya)*(y1-ya));
                rb = sqrt((x1-xb)*(x1-xb)+(y1-yb)*(y1-yb));
                ra3= ra*ra*ra;
                rb3= rb*rb*rb;
		k1[0]=dt*vx1;
		k1[1]=dt*vy1;
		k1[2]=dt*(-(x1-xa)/ra3-mb*(x1-xb)/rb3);
		k1[3]=dt*(-(y1-ya)/ra3-mb*(y1-yb)/rb3);
		x2 = x[i]+k1[0]/2.0;
		y2 = y[i]+k1[1]/2.0;
		vx2= vx[i]+k1[2]/2.0;
		vy2= vy[i]+k1[3]/2.0;
		ra = sqrt((x2-xa)*(x2-xa)+(y2-ya)*(y2-ya));
                rb = sqrt((x2-xb)*(x2-xb)+(y2-yb)*(y2-yb));
                ra3= ra*ra*ra;
                rb3= rb*rb*rb;
		k2[0]=dt*vx2;
		k2[1]=dt*vy2;
		k2[2]=dt*(-(x2-xa)/ra3-mb*(x2-xb)/rb3);
		k2[3]=dt*(-(y2-ya)/ra3-mb*(y2-yb)/rb3);
		x3 =x[i]+k2[0];
		y3 =y[i]+k2[1];
		vx3=vx[i]+k2[2];
		vy3=vy[i]+k2[3];
		ra = sqrt((x3-xa)*(x3-xa)+(y3-ya)*(y3-ya));
                rb = sqrt((x3-xb)*(x3-xb)+(y3-yb)*(y3-yb));
                ra3= ra*ra*ra;
                rb3= rb*rb*rb;
		k3[0]=dt*vx3;
		k3[1]=dt*vy3;
		k3[2]=dt*(-(x3-xa)/ra3-mb*(x3-xb)/rb3);
		k3[3]=dt*(-(y3-ya)/ra3-mb*(y3-yb)/rb3);
		x[i] =x[i] +(k0[0]+2.0*k1[0]+2.0*k2[0]+k3[0])/6.0;
                y[i] =y[i] +(k0[1]+2.0*k1[1]+2.0*k2[1]+k3[1])/6.0;
		vx[i]=vx[i]+(k0[2]+2.0*k1[2]+2.0*k2[2]+k3[2])/6.0;
		vy[i]=vy[i]+(k0[3]+2.0*k1[3]+2.0*k2[3]+k3[3])/6.0;
                }
                xa=xan;
                ya=yan;
                xb=xbn;
                yb=ybn;
                vxa=vxan;
                vya=vyan;
                vxb=vxbn;
                vyb=vybn;
        }
        printf("%f %f\n",xa,ya);
        printf("%f %f\n",xb,yb);
        for(i=0;i<N;i++)
	    printf("%f %f\n",x[i],y[i]);
}