#!/usr/bin/python3
###########################################################################################
#  package:   pNbody
#  file:      plotCylindricalProfile
#  brief:     
#  copyright: GPLv3
#             Copyright (C) 2019 EPFL (Ecole Polytechnique Federale de Lausanne)
#             LASTRO - Laboratory of Astrophysics of EPFL
#  author:    Yves Revaz <yves.revaz@epfl.ch>
#
# This file is part of Gtools.
###########################################################################################

import os
import numpy as np
import argparse
import copy

import matplotlib.pyplot as plt
from pNbody import plot
from pNbody import *

from pNbody import units, ctes

####################################################################
# option parser
####################################################################

description="""plot different physical quantities as a function of the 3-d radius, 
assuming spherical quantities."""
epilog     ="""
Examples:
--------

plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64 --log xy --center histocenter             snapshot.hdf5
plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64 --log xy --forceComovingIntegrationOff    snapshot.hdf5
plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64 --log xy --param params                   snapshot.hdf5
plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64 --log xy --remove_ic_shift                snapshot.hdf5
plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64 --log xy --select stars                   snapshot.hdf5


# use an irregular grid
plotCylindricalProfile -y sigmaz  --xmax 10 --rmax 10 --ymin 0 --nr 64 --npb1 30 --ymin 0  snap.dat

# use an irregular grid and add error bars estimated from different line of sights 
plotCylindricalProfile -y sigmaz  --xmax 10 --rmax 10 --ymin 0 --nr 64 --npb1 30 --npb2 100 --ymin 0 --nlos 50 snapshot.hdf5



# surface density profile
plotCylindricalProfile -y sdens --xmax 50 --rmax 50 --nr 64  --log xy  snapshot.hdf5

# mass profile
plotCylindricalProfile -y mass --xmax 50 --rmax 50 --nr 64  --log xy  snapshot.hdf5

# integrated mass profile
plotCylindricalProfile -y imass --xmax 50 --rmax 50 --nr 64  --log xy  snapshot.hdf5

# line of sight velocity profile
plotCylindricalProfile -y sigmaz  --xmax 50 --rmax 50 --ymin 0 --nr 64   snapshot.hdf5

# circular velocity
plotCylindricalProfile -y vcirc   --xmax 50 --rmax 50 --nr 64  --nr 64 --eps 0.1  snapshot.hdf5

"""

parser = argparse.ArgumentParser(description=description,epilog=epilog,formatter_class=argparse.RawDescriptionHelpFormatter)

plot.add_files_options(parser)
plot.add_arguments_units(parser)
plot.add_arguments_reduc(parser)
plot.add_arguments_center(parser)
plot.add_arguments_select(parser)
plot.add_arguments_info(parser)
plot.add_comoving_options(parser)
plot.add_arguments_legend(parser)
plot.add_arguments_icshift(parser)
plot.add_arguments_cmd(parser)


parser.add_argument(action="store", 
                    dest="files", 
                    metavar='FILE', 
                    type=str,
                    default=None,
                    nargs='*',
                    help='a list of files')                     


parser.add_argument("-o",
                    action="store",
                    type=str,
                    dest="outputfilename",
                    default=None,
                    help="Name of the output file")  



parser.add_argument('--xmin',
                    action="store", 
                    dest="xmin", 
                    metavar='FLOAT', 
                    type=float,
                    default=None,
                    help='x min')

parser.add_argument('--xmax',
                    action="store", 
                    dest="xmax", 
                    metavar='FLOAT', 
                    type=float,
                    default=None,
                    help='x max')
                    
parser.add_argument('--ymin',
                    action="store", 
                    dest="ymin", 
                    metavar='FLOAT', 
                    type=float,
                    default=None,
                    help='y min')

parser.add_argument('--ymax',
                    action="store", 
                    dest="ymax", 
                    metavar='FLOAT', 
                    type=float,
                    default=None,
                    help='y max')

                                        
parser.add_argument('--log',
                    action="store", 
                    dest="log", 
                    metavar='STR', 
                    type=str,
                    default=None,
                    help='log scale (None,x,y,xy)')


parser.add_argument('-y','--y',
                    action="store", 
                    dest="y", 
                    metavar='STR', 
                    type=str,
                    default='density',
                    help='quantity to plot in the y axis')

parser.add_argument("--rmax",
                  action="store",
                  dest="rmax",
                  type=float,
                  default=50.,
                  help="max radius of bins",
                  metavar=" FLOAT")

parser.add_argument("--nr",
                  action="store",
                  dest="nr",
                  type=int,
                  default=32,
                  help="number of bins in r",
                  metavar=" INT")   

parser.add_argument("--nt",
                     action="store",
                     dest="nt",
                     type=int,
                     default=32,
                     help="number of division in theta",
                     metavar=" INT")
                     
parser.add_argument("--eps",
                     action="store",
                     dest="eps",
                     type=float,
                     default=0.1,
                     help="smoothing length",
                     metavar=" FLOAT")
   
   
parser.add_argument("--npb1",
                     action="store",
                     dest="npb1",
                     type=int,
                     default=0,
                     help="number of particles per bin (sub plots)",
                     metavar=" INT")
   
parser.add_argument("--npb2",
                     action="store",
                     dest="npb2",
                     type=int,
                     default=30,
                     help="number of particles per bin (total plot)",
                     metavar=" INT")
   

parser.add_argument("--nlos",
                      action="store",
                      dest="nlos",
                      type=int,
                      default=0,
                      help="number of line of sight to use",
                      metavar=" INT")

parser.add_argument("--rmode",
                      action="store",
                      dest="rmode",
                      type=str,
                      default='uniform',
                      help="mode used to define the bin radius in irregular grid (center, mean, uniform)",
                      metavar=" INT")


parser.add_argument("--AdaptativeSoftenning",
                      action="store_true",
                      dest="AdaptativeSoftenning",
                      default=False,
                      help="AdaptativeSoftenning")

parser.add_argument("--ErrTolTheta",
                      action="store",
                      dest="ErrTolTheta",
                      type=float,
                      default=0.5,
                      help="Error tolerance theta",
                      metavar=" FLOAT")



parser.add_argument("--colormap",
                    action="store", 
                    default='jet',
                    help='matplotlib colormap name (e.g. mycmap, tab20c, Greys, jet, binary)') 
                    
             


#######################################
# MakePlot
#######################################


def MakePlot(opt):
  
  params = {
    "axes.labelsize": 14,
    "axes.titlesize": 18,
    "font.size": 12,
    "legend.fontsize": 12,
    "xtick.labelsize": 14,
    "ytick.labelsize": 14,
    "text.usetex": True,
    "figure.subplot.left": 0.15,
    "figure.subplot.right": 0.95,
    "figure.subplot.bottom": 0.15,
    "figure.subplot.top": 0.95,
    "figure.subplot.wspace": 0.02,
    "figure.subplot.hspace": 0.02,
    "figure.figsize" : (8, 6),
    "lines.markersize": 6,
    "lines.linewidth": 2.0,
  }
  plt.rcParams.update(params)
  
  
  # create the plot
  fig = plt.gcf()
  fig.set_size_inches(8,6)
  ax  = plt.gca()


  # get a list of color
  colors = plot.ColorList(n=len(opt.files),colormap=opt.colormap)
  
  # list of points
  datas = []
  
  
  #################################
  # loop over files
  
  for filename in opt.files:
    
    nb = Nbody(filename, ftype=opt.ftype)
    
    ################
    # units
    ################
    
    # define local units
    unit_params = plot.apply_arguments_units(opt)
    nb.set_local_system_of_units(params=unit_params)
    
    ################
    # apply options
    ################
    nb = plot.apply_arguments_icshift(nb, opt)
    nb = plot.apply_arguments_comoving(nb, opt)
    nb = plot.apply_arguments_reduc(nb, opt)
    nb = plot.apply_arguments_select(nb, opt)
    nb = plot.apply_arguments_center(nb, opt)
    nb = plot.apply_arguments_cmd(nb, opt)
    nb = plot.apply_arguments_info(nb, opt)
    nb = plot.apply_arguments_display(nb, opt)


    ################
    # some info
    ################
    print("---------------------------------------------------------")
    nb.localsystem_of_units.info()
    nb.ComovingIntegrationInfo()
    print("---------------------------------------------------------")

      

    # grid division
    rc = 1

    def f(r): return np.log(r / rc + 1.)

    def fm(r): return rc * (np.exp(r) - 1.)

    ###############################
    # compute physical quantities
    ###############################
    
    
    if opt.npb1 == 0:

      ##########################################
      # use Cylindrical_2drt_Grid
      ##########################################
      
      if opt.y == 'sdens':
      
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=1, g=f, gm=fm)
      
          x, t = G.get_rt()
          y = G.get_SurfaceDensityMap(nb)
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
              x = x * nb.atime / nb.hubbleparam		# length  conversion
              y = y / nb.atime**2 * nb.hubbleparam 	# surface density conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitSurfaceDensity)
      
          x = x * fx
          y = y * fy
      
          # set labels
      
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Surface\,\,Density}\,\left[ M_{\odot}/kpc^2 \right]$'
      
          x, y = plot.CleanVectorsForLogX(x, y)
          x, y = plot.CleanVectorsForLogY(x, y)
          datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      
      if opt.y == 'mass':
      
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=1, g=f, gm=fm)
      
          x, t = G.get_rt()
          y = G.get_MassMap(nb)
          #y = sum(y, axis=1)
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
              x = x * nb.atime / nb.hubbleparam		# length  conversion
              y = y / nb.hubbleparam 			# mass conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitDensity)
      
          x = x * fx
          y = y * fy
      
          # set labels
      
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$'
      
          x, y = plot.CleanVectorsForLogX(x, y)
          x, y = plot.CleanVectorsForLogY(x, y)
          datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      
      if opt.y == 'imass':
      
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=1, g=f, gm=fm)
      
          x, t = G.get_rt()
          y = G.get_MassMap(nb)
          #y = sum(y, axis=1)
          y = np.add.accumulate(y)
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
              x = x * nb.atime / nb.hubbleparam		# length  conversion
              y = y / nb.hubbleparam 			# mass conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitDensity)
      
          x = x * fx
          y = y * fy
      
          # set labels
      
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Mass}\,\left[ M_{\odot} \right]$'
      
          x, y = plot.CleanVectorsForLogX(x, y)
          x, y = plot.CleanVectorsForLogY(x, y)
          datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      
      if ((opt.y == 'sigmar') or (opt.y == "sigmat") or (opt.y == "sigmaz")) :
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
      
              x0 = nb.pos
              v0 = nb.vel
      
              Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
              pars = {
                  "Hubble": Hubble,
                  "HubbleParam": nb.hubbleparam,
                  "OmegaLambda": nb.omegalambda,
                  "Omega0": nb.omega0}
              a = nb.atime
              Ha = cosmo.Hubble_a(a, pars=pars)
      
              nb.pos = x0 * nb.atime / nb.hubbleparam  # length    conversion
              nb.vel = v0 * np.sqrt(a) + x0 * Ha * a		# velocity  conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_km, units.Unit_Ms, units.Unit_s, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitVelocity)
      
          nb.pos = nb.pos * fx
          nb.vel = nb.vel * fy
          
          
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=1, g=f, gm=fm)
      
          x, t = G.get_rt()
          
          if opt.y == 'sigmaz':
            y = G.get_SigmaValMap(nb, nb.Vz())
            ylabel = r'$\sigma_z\,\left[ \rm{km}/\rm{s} \right]$'
      
          if opt.y == 'sigmar':
            y = G.get_SigmaValMap(nb, nb.VR())
            ylabel = r'$\sigma_R\,\left[ \rm{km}/\rm{s} \right]$'              
      
          if opt.y == 'sigmat':
            y = G.get_SigmaValMap(nb, nb.VT())
            ylabel = r'$\sigma_\theta\,\left[ \rm{km}/\rm{s} \right]$' 
      

      
          # set labels
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
      
      
      
          
          if opt.nlos > 1:
            G = libgrid.CylindricalIrregular_1dr_Grid(xs,opt.npb2)
            x,mean,std = G.get_ValMeanAndStd(ys)
            y = mean
            
            datas.append(pt.DataPoints(x,y,color=colors.get(),label=filename,tpe='points'))
            datas.append(pt.DataPoints(x,y+std,color='r',label=filename,tpe='points'))
            datas.append(pt.DataPoints(x,y-std,color='r',label=filename,tpe='points'))
            
            
      
          else:
            
            x, y = plot.CleanVectorsForLogX(x, y)
            x, y = plot.CleanVectorsForLogY(x, y)
                        
            datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      
      
      
      
      if opt.y == 'vcirc':
      
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=opt.nt, z=0, g=f, gm=fm)
      
          x, t = G.get_rt()
      
          # compute tree
          nb.getTree(force_computation=True, ErrTolTheta=opt.ErrTolTheta)
          # compute accel
          Accx, Accy, Accz = G.get_AccelerationMap(
              nb, eps=opt.eps, UseTree=True, AdaptativeSoftenning=opt.AdaptativeSoftenning)
      
          Ar = np.sqrt(Accx**2 + Accy**2)
          Ar = np.sum(Ar, 1) / opt.nt
          dPhit = Ar
          
          
          y = libdisk.Vcirc(x, dPhit)
          # correct from Gravity constant
          Gcte = ctes.GRAVITY.into(nb.localsystem_of_units)
          y = y*np.sqrt(Gcte)
          
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
              x = x * nb.atime / nb.hubbleparam		# length  conversion
              # y = y/nb.hubbleparam 			# velocity
              print("!!! this is not defined for comobiles coords, please, do it !!!")
              sys.exit()
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_km, units.Unit_Ms, units.Unit_s, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitVelocity)
      
          x = x * fx
          y = y * fy
      
          # set labels
      
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Circular\,Velocity}\,\left[ \rm{km}/\rm{s} \right]$'
      
          x, y = plot.CleanVectorsForLogX(x, y)
          x, y = plot.CleanVectorsForLogY(x, y)
          datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      
      if opt.y == 'Lum':
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
      
              x0 = nb.pos
              v0 = nb.vel
      
              Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
              pars = {
                  "Hubble": Hubble,
                  "HubbleParam": nb.hubbleparam,
                  "OmegaLambda": nb.omegalambda,
                  "Omega0": nb.omega0}
              a = nb.atime
              Ha = cosmo.Hubble_a(a, pars=pars)
      
              nb.pos = x0 * nb.atime / nb.hubbleparam  # length    conversion
              nb.vel = v0 * np.sqrt(a) + x0 * Ha * a		# velocity  conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_km, units.Unit_Ms, units.Unit_s, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitVelocity)
      
          nb.pos = nb.pos * fx
          nb.vel = nb.vel * fy
      
          G = libgrid.Cylindrical_2drt_Grid(
              rmin=0, rmax=opt.rmax, nr=opt.nr, nt=1, g=f, gm=fm)
      
          x, t = G.get_rt()
          y = np.reshape(
              G.get_SurfaceValMap(
                  nb,
                  nb.LuminositySpec()),
              opt.nr)		# in Lsol/(unitlength)
      
          print("the surface density units is wrong here...!!!")
      
          # set labels
      
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Luminosity}\left[ L_{\odot}/\rm{pc}^2 \right]$'
      
          if opt.log:
              x, y = plot.CleanVectorsForLogX(x, y)
              x, y = plot.CleanVectorsForLogY(x, y)
      
          datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))



    else:

      ##########################################
      # use CylindricalIrregular_1dr_Grid
      ##########################################   


      if opt.y == 'sdens':
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
              x = x * nb.atime / nb.hubbleparam		# length  conversion
              y = y / nb.atime**2 * nb.hubbleparam 	# surface density conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitSurfaceDensity)
      
      
      
          xs = np.zeros(0)
          ys = np.zeros(0)
          
          if opt.nlos>1:      
            los = plot.getLOS(opt.nlos,seed=0)
          else:
            los = np.array([0,0,1])
            
          for j in range(len(los)):
                            
            # copy the model
            nbc = copy.deepcopy(nb)
      
            # rotate
            if opt.nlos > 1:
              nbc.align(axis=los[j])
              
            # compute the radius
            r0 = nbc.rxy()
        
            # define the grid
            G = libgrid.CylindricalIrregular_1dr_Grid(r0,opt.npb1,rmode=opt.rmode)
            x = G.get_R()
            y = G.get_SurfaceDensity(nbc.mass)
            
            x = x * fx
            y = y * fy
                
            xs     = np.concatenate((xs,x))
            ys     = np.concatenate((ys,y))          
      
      

          if opt.nlos > 1:
            G = libgrid.CylindricalIrregular_1dr_Grid(xs,opt.npb2)
            x = G.get_R()
            mean,std = G.get_MeanAndStd(ys)
            y = mean
            
            datas.append(plot.DataPoints(x,y,color=colors.get(),label=filename,tpe='points'))
            datas.append(plot.DataPoints(x,y+std,color='r',label=filename,tpe='points'))
            datas.append(plot.DataPoints(x,y-std,color='r',label=filename,tpe='points'))

          else:
            x, y = plot.CleanVectorsForLogX(x, y)
            x, y = plot.CleanVectorsForLogY(x, y)
            datas.append(plot.DataPoints(x,y,color=colors.get(),label=filename,tpe='points'))
          

      
          # set labels
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
          ylabel = r'$\rm{Surface\,\,Density}\,\left[ M_{\odot}/kpc^2 \right]$'
          
          
      







      if ((opt.y == 'sigmar') or (opt.y == "sigmat") or (opt.y == "sigmaz")) :
      
          # comoving conversion
          if nb.isComovingIntegrationOn():
              print((
                  "    converting to physical units (a=%5.3f h=%5.3f)" %
                  (nb.atime, nb.hubbleparam)))
      
              x0 = nb.pos
              v0 = nb.vel
      
              Hubble = ctes.HUBBLE.into(nb.localsystem_of_units)
              pars = {
                  "Hubble": Hubble,
                  "HubbleParam": nb.hubbleparam,
                  "OmegaLambda": nb.omegalambda,
                  "Omega0": nb.omega0}
              a = nb.atime
              Ha = cosmo.Hubble_a(a, pars=pars)
      
              nb.pos = x0 * nb.atime / nb.hubbleparam  # length    conversion
              nb.vel = v0 * np.sqrt(a) + x0 * Ha * a		# velocity  conversion
      
          # output units
          out_units_x = units.UnitSystem(
              'local', [units.Unit_kpc, units.Unit_Ms, units.Unit_Myr, units.Unit_K])
          out_units_y = units.UnitSystem(
              'local', [units.Unit_km, units.Unit_Ms, units.Unit_s, units.Unit_K])
      
          fx = nb.localsystem_of_units.convertionFactorTo(
              out_units_x.UnitLength)
          fy = nb.localsystem_of_units.convertionFactorTo(
              out_units_y.UnitVelocity)
      
          nb.pos = nb.pos * fx
          nb.vel = nb.vel * fy
          
            
          xs = np.zeros(0)
          ys = np.zeros(0)
          
          if opt.nlos>1:      
            los = plot.getLOS(opt.nlos,seed=0)
          else:
            los = np.array([0,0,1])
            
          for j in range(len(los)):
                            
            # copy the model
            nbc = copy.deepcopy(nb)
      
            # rotate
            if opt.nlos > 1:
              nbc.align(axis=los[j])
              
            # compute the radius
            r0 = nbc.rxy()
        
            # compute the variable 
            if opt.y == 'sigmaz':
              val = nbc.Vz()
            if opt.y == 'sigmar':
              val = nbc.VR()
            if opt.y == 'sigmat':
              val= nbc.VT()
            
            # define the grid
            G = libgrid.CylindricalIrregular_1dr_Grid(r0,opt.npb1,rmode=opt.rmode)
      
            # compute the values
            rs = G.get_R()
            Means,Sigmas = G.get_MeanAndStd(val)
            
            x = rs 
            y = Sigmas
      
            xs     = np.concatenate((xs,x))
            ys     = np.concatenate((ys,y))
      
      
      
      
          if opt.y == 'sigmaz':
            ylabel = r'$\sigma_z\,\left[ \rm{km}/\rm{s} \right]$'
      
          if opt.y == 'sigmar':
            ylabel = r'$\sigma_R\,\left[ \rm{km}/\rm{s} \right]$'              
      
          if opt.y == 'sigmat':
            ylabel = r'$\sigma_\theta\,\left[ \rm{km}/\rm{s} \right]$' 
              
      
          # set labels
          xlabel = r'$\rm{Radius}\,\left[ \rm{kpc} \right]$'
      

          if opt.nlos > 1:
            G = libgrid.CylindricalIrregular_1dr_Grid(xs,opt.npb2)
            x = G.get_R()
            mean,std = G.get_MeanAndStd(ys)
            y = mean
            
            datas.append(plot.DataPoints(x,y,color=colors.get(),label=filename,tpe='points'))
            datas.append(plot.DataPoints(x,y+std,color='r',label=filename,tpe='points'))
            datas.append(plot.DataPoints(x,y-std,color='r',label=filename,tpe='points'))
            

          else:
            
            x, y = plot.CleanVectorsForLogX(x, y)
            x, y = plot.CleanVectorsForLogY(x, y)
                        
            datas.append(
              plot.DataPoints(
                  x,
                  y,
                  color=colors.get(),
                  label=filename,
                  tpe='points'))
      


          
  ##################
  # plot all
  ##################

  for d in datas:
    ax.plot(d.x, d.y, color=d.color)

  # set limits
  xmin, xmax, ymin, ymax, log = plot.SetLimitsFromDataPoints(datas, opt.xmin, opt.xmax, opt.ymin, opt.ymax, opt.log)

  # set the axis (the extention is done in the previous command)
  plot.SetAxis(ax,xmin,xmax,ymin,ymax,log,extend=False)
  
  # labels
  ax.set_xlabel(xlabel)
  ax.set_ylabel(ylabel)
  
  # legend
  if opt.legend:
    plot.LegendFromDataPoints(ax, datas, opt.legend_loc)



  # save or display
  if opt.outputfilename:
    plt.savefig(opt.outputfilename)
  else:
    plt.show()    
      


#################################
# main
#################################

if __name__ == '__main__':  
  
  opt = parser.parse_args()
  MakePlot(opt)
