  We propose a new numerical method to calculate irrotational binary systems
composed of compressible gaseous stars in Newtonian gravity. Assuming
irrotationality, i.e. vanishing of the vorticity vector everywhere in the star
in the inertial frame, we can introduce the velocity potential for the flow
field. Using this velocity potential we can derive a set of basic equations for
stationary states which consist of (i) the generalized Bernoulli equation, (ii)
the Poisson equation for the Newtonian gravitational potential and (iii) the
equation for the velocity potential with the Neumann type boundary condition.
We succeeded in developing a new code to compute numerically exact solutions to
these equations for the first time. Such irrotational configurations of binary
systems are appropriate models for realistic neutron star binaries composed of
inviscid gases, just prior to coalescence of two stars caused by emission of
gravitational waves. Accuracies of our numerical solutions are so high that we
can compute reliable models for fully deformed final stationary configurations
and hence determine the inner most stable circular orbit of binary neutron star
systems under the approximations of weak gravity and inviscid limit.
