Planetary orbits

c planetary orbits using Cartesian coordinates
      real M   ,mearth
      data  M ,mearth,G , x0,y0, v0/ 1.98E30,5.98e24, 6.67e-11,
     $ 1.50E11,0.,2.99E4/
      ax(x,y)=-G*m*x/(x**2+y**2)**(3./2.)
      ay(x,y)=-G*m*y/(x**2+y**2)**(3./2.)
      pi=2.*asin(1.)
      alfa=90.
      alfa=alfa*pi/180.
      tscale=sqrt(x0**2+y0**2)/v0
      dt=tscale/100.
      tfinal=6.*tscale
      nstep=int(tfinal/dt)
      kp=int(float(nstep)/60.)
      kount=kp
      v0x=v0*cos(alfa)
      v0y=v0*sin(alfa)
      x1=x0+v0x*dt
      y1=y0+v0y*dt
      Ekin=.5*mearth*(v0x**2+v0y**2)
      Epot=-G*m*mearth/(x0**2+y0**2)**(.5)
      etotal=Ekin+epot
      print*,'etotal(j)=', etotal
      Print*,'   '
      print 100,0., x0, y0
      do 10 i=2,nstep
      t=dt*float(i)
      x2=2.*x1-x0+dt**2*ax(x1,y1)
      y2=2.*y1-y0+dt**2*ay(x1,y1)
      vx=(x2-x1)/dt
      vy=(y2-y1)/dt
      if(i.eq.kount)then
      ekin=.5*mearth*(vx**2+vy**2)
      Epot=-G*m*mearth/(x2**2+y2**2)**(.5)
      etotal=ekin+epot
      print 100,t, x2, y2
      kount=kount+kp
      endif
      x0=x1
      x1=x2
      y0=y1
      y1=y2
10    continue
100   format(1x,'t,x,y=',3(4x,e10.3) )
      stop
      end