function corte_lin,teta0,b,ns,xa,ya,za,x0s,del0,dirz,z0,xb,yb,zb,p0
  Nl=fltarr(ns)
  Nlp=fltarr(ns)
  Nln=fltarr(ns)
  m=tan(teta0*!pi/180.)
  xl=fltarr(ns)
  yl=fltarr(ns)
  rl=fltarr(ns)
  rl0=fltarr(ns)
  Bxl=fltarr(ns)
  Byl=fltarr(ns)
  Bzl=fltarr(ns)
  BBl=fltarr(ns)
  
  print,' working with the angle = ',teta0

  x0=x0s[0]
  y0=x0*m+b
  for i=0,ns-1 do begin
     xi=x0s[i]
     yi=xi*m+b
     xl[i]=xi
     yl[i]=yi     ;aplicando la recta y calculando la dist al punto mas cercano
     dxi=xi-p0[0]
     dyi=yi-p0[1]
     rl[i]=sqrt((dxi)^2+(dyi)^2); la dist 
    ; indb=where(sqrt((xa - xi)^2 + (ya - yi)^2 + (za)^2) le del0/2.,countb)
     indb=where(sqrt((xa - xi)^2 + (ya - yi)^2) le del0/10.,countb) ;para el campo se divide entre 10
     if countb ge 1 then begin 
  ;      print,i,countb
        Bxl[i]=mean(xb[indb])
        Byl[i]=mean(yb[indb])
        Bzl[i]=mean(zb[indb])
        BBl[i]=sqrt(Bxl[i]^2+ Byl[i]^2+Bzl[i]^2)
     endif
     ; print,'dist ',x0,xi,y0,yi,rl[i]
     nae=0
     naep=0
     naen=0
     sz=size(xa)   ;el numero de trayectorias 
     for k=0, sz[1]-1 do begin ;iteramos sobre cada trayctoria por separado
        indi=where( sqrt((xa[k,*] - xi)^2 + (ya[k,*] - yi)^2 + (za[k,*] - Z0)^2) le del0,countk)
         if countk  gt 0 then begin
            tdz=total((dirz[k,indi]))     ;    eq 1,countip)   
            nae  = nae  + 1               ;counti ;sumamos uno cuado toca la
            if tdz gt 0 then naep = naep + 1 ; countip ;esfera !!!!
            if tdz lt 0 then naen = naen + 1 ; countin
 ;           print,' k = ',k,countk
         endif
     endfor
     Nl[i]=nae    
     Nlp[i]=naep
     Nln[i]= naen
  endfor

  disti=rl                      ;sqrt(dxi^2+dyi^2)
  dind=where(disti eq min(disti))
  rl=[(disti[0:dind[0]-1]-disti[dind[0]])*(-1),disti[dind[0]:-1]-disti[dind[0]]]
  sal=[[rl],[nl],[nlp],[nln],[xl],[yl],[Bxl],[Byl],[Bzl],[BBl]]
 return, sal
end

pro plot_radial_profile15_gen,dir,name,se,dang,distancia,graf,bufer
 
;dang in this case maybe -50, 0 or 50
;  Vamos a calcular una esfera del diametro de la Tierra
  CPU, TPOOL_MIN_ELTS=5000, TPOOL_NTHREADS=3
   if graf ne 0 then loadct,12
  if bufer eq 1 then buf=1 else buf=0 
                                ; basico_sim_fr
restore,'/home/alara/ARTS/ANDREA/SIM/DATA/basico_sim.save'
sname='data_'+name+'_all.save'  ;'data_eduardo_orig.save'

p0s=fltarr(2,3)
sse=strcompress(se,/rem)
sdist=strcompress(distancia,/rem)
  case dang of
     -50 : begin
        tetals=tetals50N
        fdir=dir+'50N'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_50N.save'
        signo=1
     end
     0   : begin
        tetals=tetals0
        fdir=dir+'0'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_0.save'
        signo=1
        end
     50 : begin
        tetals=tetals50
        fdir=dir+'50'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_50.save'
        signo=-1
      end
     -20 : begin
        tetals=tetals20N
        fdir=dir+'20N'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_20N.save'
        signo=1
     end
     20 : begin
        tetals=tetals20
        fdir=dir+'20'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_20.save'
        signo=-1
        end
     -15 : begin
        tetals=tetals15N
        fdir=dir+'15N'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_15N.save'
        signo=1
     end
     15 : begin
        tetals=tetals15
        fdir=dir+'15'
        ssname='sim_'+name+'_'+sse+'_'+sdist+'_15.save'
        signo=-1
        end
   endcase

  help,distancia,tetals
  print,distancia,tetals
  
        p0s[*,*]=[[distancia*cos((90-tetals)*!pi/180.)],[distancia*sin((90-tetals)*!pi/180.)]]
        b0s=signo*abs(distancia/sin((90-tetals)*!pi/180.))
        print,'b0s = ',b0s
        print,'p0s = ',p0s

  
    
ddir='/home/alara/ARTS/ANDREA/SIM/DATA/'+dir+'/' ;OCT16/dest_folder/'
fdir='/home/alara/ARTS/ANDREA/SIM/FIG/'+fdir+'/'  ;OCT16/dest_folder/'


  n0=0.1/(10^(-2.7*11.))

  fils=file_search(ddir,sname)
  print,'file = ',fils
  restore,fils
  ;help,lsx,lsy,lsz,lse
  ax=LSX.ToArray()
  ay=LSy.ToArray()
  az=LSZ.ToArray()
  ae=LSE.ToArray()
  aalpha=LSALPHA.ToArray()
  at=LST.ToArray()
  adirz=LSdirZ.ToArray()
   if graf ne 0 then colores= !color
  bx=LSBX.ToArray()
  by=LSBy.ToArray()
  bz=LSBz.ToArray()

 sae=sort(ae)
 indeu=uniq(ae[sae])
 es=ae[sae[indeu]] 
 nne=n_elements(es)
 nes=es^(-2.7)*10d^2.7

  ;para la division en z
  nsz=  25
  zmin=-0.25
  zmax=0.25
  z0s=indgen(nsz)*(zmax - zmin)/(nsz-1) +zmin

                                
  linl20=list()
  linl0=list()
  linl15=list()
  linl30=list()
  linl20p=list()
  linl0p=list()
  linl15p=list()
  linl30p=list()
  linl20n=list()
  linl0n=list()
  linl15n=list()
  linl30n=list()
  linl55=list()
  linl55p=list()
  linl55n=list()

       ;cambiando alturas
  for j=0,nsz-1 do begin
     ; zri=az+z0s[j]
     print,'working with z = ',z0s[j]
     inde=where(ae eq se)
     print,'B0s = ',b0s[0],b0s[1],b0s[2]
     print,'P0s = ',p0s
     lini20=corte_lin(tetals[0],b0s[0],ns,ax[inde,*],ay[inde,*],az[inde,*],x0s,del0,adirz[inde,*],z0s[j],bx[inde,*],by[inde,*],bz[inde,*],p0s[*,0])
     ;print,'lin20 = ',lini20
     dsal20=reform(lini20[*,0])
     linl20.Add,reform(lini20[*,1])
     linl20p.Add,reform(lini20[*,2])
     linl20n.Add,reform(lini20[*,3])
     xsal20=reform(lini20[*,4])  ; suponemos que no cambian x y y en la iter 
     ysal20=reform(lini20[*,5])
     Bxsal20=reform(lini20[*,6])   
     Bysal20=reform(lini20[*,7])
     Bzsal20=reform(lini20[*,8])   
     BBsal20=reform(lini20[*,9])
    
     lini0=corte_lin(tetals[1],b0s[1],ns,ax[inde,*],ay[inde,*],az[inde,*],x0s,del0,adirz[inde,*],z0s[j],bx[inde,*],by[inde,*],bz[inde,*],p0s[*,1])
     dsal0=reform(lini0[*,0])
     linl0.Add,reform(lini0[*,1])
     linl0p.Add,reform(lini0[*,2])
     linl0n.Add,reform(lini0[*,3])
     xsal0=reform(lini0[*,4])
     Bxsal0=reform(lini0[*,6])   
     Bysal0=reform(lini0[*,7])
     Bzsal0=reform(lini0[*,8])   
     BBsal0=reform(lini0[*,9])
     ysal0=reform(lini0[*,5])

     lini15=corte_lin(tetals[2],b0s[2],ns,ax[inde,*],ay[inde,*],az[inde,*],x0s,del0,adirz[inde,*],z0s[j],bx[inde,*],by[inde,*],bz[inde,*],p0s[*,2])
     dsal15=reform(lini15[*,0])
     linl15.Add,reform(lini15[*,1])
     linl15p.Add,reform(lini15[*,2])
     linl15n.Add,reform(lini15[*,3])
     xsal15=reform(lini15[*,4])
     ysal15=reform(lini15[*,5])
     Bxsal15=reform(lini15[*,6])   
     Bysal15=reform(lini15[*,7])
     Bzsal15=reform(lini15[*,8])   
     BBsal15=reform(lini15[*,9])
 
endfor

  if graf ne 0 then begin
  pltt=plot(xsal20,ysal20,buffer=buf)
  pltt=plot(xsal0,ysal0,'b',/over)
  pltt=plot(xsal15,ysal15,'g',/over)
  pltt.Close
  endif
  
la20=total(linl20.ToArray(),1)
la0=total(linl0.ToArray(),1)
la15=total(linl15.ToArray(),1)
;la30=total(linl30.ToArray(),1)
;la55=total(linl55.ToArray(),1)

la20p=total(linl20p.ToArray(),1)
la0p=total(linl0p.ToArray(),1)
la15p=total(linl15p.ToArray(),1)
;la30p=total(linl30p.ToArray(),1)
;la55p=total(linl55p.ToArray(),1)

la20n=total(linl20n.ToArray(),1)
la0n=total(linl0n.ToArray(),1)
la15n=total(linl15n.ToArray(),1)
;la30n=total(linl30n.ToArray(),1)
;la55n=total(linl55n.ToArray(),1)

ind20ri=where(dsal20 lt -riny or dsal20 gt riny,count20ri)
ind0ri=where(dsal0 lt -riny or dsal0 gt riny,count0ri)
 ind15ri=where(dsal15 lt -riny or dsal15 gt riny,count15ri)
;ind30ri=where(dsal30 lt -riny or dsal30 gt riny,count30ri)
;nout=total(nes)


save,filename=ddir+ssname,la20,la0,la15,dsal20,dsal0,dsal15,la20p,la0p,la15p,la20n,la0n,la15n,se,xsal20,ysal20,xsal0,ysal0,xsal15,ysal15,Bxsal20,Bysal20,Bzsal20,BBsal20,Bxsal0,BBsal0,Bysal0,Bzsal0,Bxsal15,Bysal15,Bzsal15,BBsal15
 print,ssname

 
  if graf ne 0 then begin
figname=fdir+'sim_'+name+'_'+sse+'_'+sdist+'.png'
plt1=plot(dsal20,smooth(LA20,15),'g',xtitl='distance (AU)',ytitl='N',titl='E = '+strcompress(se,/rem),buffer=buf);,xrange=[-0.4,0.4])
lt1=plot(dsal0,smooth(LA0,15),/over)
plt1=plot(dsal15,smooth(LA15,15),'b',/over)
;plt1=plot(dsal30,smooth(LA30,15),'magenta',/over)
;plt1=plot(dsal55,smooth(LA55,15),'orange',/over)
t1=text(0.55,0.84,strcompress(tetals[0]),'g')
t1=text(0.55,0.8,strcompress(tetals[1]),'black')
t1=text(0.55,0.76,strcompress(tetals[2]),'b')
;t1=text(0.55,0.72,strcompress(tetals[3]),'magenta')
;t1=text(0.55,0.68,strcompress(tetals[4]),'orange')
plt1.Save,figname,res=300
plt1.Close

figname2=fdir+'sim_p_'+name+'_'+sse+'_'+sdist+'.png'
plt2=plot(dsal20,smooth(LA20p,15),'g',xtitl='distance (AU)',ytitl='N',titl='E = '+strcompress(se,/rem)+' +z',buffer=buf);,xrange=[-0.4,0.4])
lt1=plot(dsal0,smooth(LA0p,15),/over)
plt2=plot(dsal15,smooth(LA15p,15),'b',/over)
;plt2=plot(dsal30,smooth(LA30p,15),'magenta',/over)
;plt2=plot(dsal55,smooth(LA55p,15),'orange',/over)
t1=text(0.55,0.84,strcompress(tetals[0]),'g')
t1=text(0.55,0.8,strcompress(tetals[1]),'black')
t1=text(0.55,0.76,strcompress(tetals[2]),'b')
;t1=text(0.55,0.72,strcompress(tetals[3]),'magenta')
;t1=text(0.55,0.68,strcompress(tetals[4]),'orange')
plt2.Save,figname2,res=300
plt2.Close

figname3=fdir+'sim_n_'+name+'_'+sse+'_'+sdist+'.png'
plt3=plot(dsal20,smooth(LA20n,15),'g',xtitl='distance (AU)',ytitl='N',titl='E = '+strcompress(se,/rem)+' -z',buffer=buf);,xrange=[-0.4,0.4])
lt1=plot(dsal0,smooth(LA0n,15),/over)
plt3=plot(dsal15,smooth(LA15n,15),'b',/over)
;plt3=plot(dsal30,smooth(LA30n,15),'magenta',/over)
;plt3=plot(dsal55,smooth(LA55n,15),'orange',/over)
t1=text(0.55,0.84,strcompress(tetals[0]),'g')
t1=text(0.55,0.8,strcompress(tetals[1]),'black')
t1=text(0.55,0.76,strcompress(tetals[2]),'b')
;t1=text(0.55,0.72,strcompress(tetals[3]),'magenta')
;t1=text(0.55,0.68,strcompress(tetals[4]),'orange')
plt3.Save,figname3,res=300

plt3.Close
endif
end