subroutine script_common(scriptunit, areab1b2, square, seedname)
!================================================!
use w90_constants, only: dp
integer, intent(in) :: scriptunit
real(kind=dp), intent(in) :: areab1b2
character(len=25), intent(in) :: square
character(len=50), intent(in) :: seedname
write (scriptunit, '(a)') "import pylab as pl"
write (scriptunit, '(a)') "import numpy as np"
write (scriptunit, '(a)') "from scipy import interpolate"
write (scriptunit, '(a)') "from collections import OrderedDict"
write (scriptunit, '(a)') " "
write (scriptunit, '(a)') "points = np.loadtxt('"//trim(seedname)// &
"-kslice-coord.dat')"
write (scriptunit, '(a)') "# Avoid numerical noise"
write (scriptunit, '(a)') "points_x=np.around(points[:,0],decimals=10)"
write (scriptunit, '(a)') "points_y=np.around(points[:,1],decimals=10)"
write (scriptunit, '(a)') "num_pt=len(points)"
write (scriptunit, '(a)') " "
write (scriptunit, '(a,f12.6)') "area=", areab1b2
write (scriptunit, '(a)') " "
write (scriptunit, '(a)') "square= "//square
write (scriptunit, '(a)') " "
write (scriptunit, '(a)') "if square:"
write (scriptunit, '(a)') &
" x_coord=list(OrderedDict.fromkeys(points_x))"
write (scriptunit, '(a)') &
" y_coord=list(OrderedDict.fromkeys(points_y))"
write (scriptunit, '(a)') " dimx=len(x_coord)"
write (scriptunit, '(a)') " dimy=len(y_coord)"
write (scriptunit, '(a)') "else:"
write (scriptunit, '(a)') " xmin=np.min(points_x)"
write (scriptunit, '(a)') " ymin=np.min(points_y)"
write (scriptunit, '(a)') " xmax=np.max(points_x)"
write (scriptunit, '(a)') " ymax=np.max(points_y)"
write (scriptunit, '(a)') &
" a=np.max(np.array([xmax-xmin,ymax-ymin]))"
write (scriptunit, '(a)') &
" num_int=int(round(np.sqrt(num_pt*a**2/area)))"
write (scriptunit, '(a)') " xint = np.linspace(xmin,xmin+a,num_int)"
write (scriptunit, '(a)') " yint = np.linspace(ymin,ymin+a,num_int)"
write (scriptunit, '(a)') " "
end subroutine script_common