;# cora001, cora002

pi = acos(-1.)
nz = 501.
zmax = 3.
z1 = zmax/(nz-1)
zmid = zmax/2.
z = z1*findgen(nz)-zmid

zmx21 = 1.8
zmx22 = 3.0
edge  = 0.15
tzr   = 0.2 
freq  = 1.
;# seycos
  arg = 2.*pi/zmax 
  norm = zmax/(2.*pi)
  zpi = 2.*pi/zmax 
  zinv = 1./zpi

  t2 = (.5*freq+tzr*tzr)*z+.5*freq*(zinv*sin(arg*z)+ 1.)                
  n2 = (.5*freq+tzr*tzr)  +.5*freq*(     cos(arg*z)    )                

;# seyduct          
  da1 = fltarr(nz)
  da2 = fltarr(nz)
  da3 = fltarr(nz)
  da4 = fltarr(nz)

  a1 = fltarr(nz)
  a2 = fltarr(nz)
  a3 = fltarr(nz)
  a4 = fltarr(nz)

  zmx1 = zmax-zmx21
  zc = z - zmx21
  zc = z - zmx21/2.
  f1 = (tzr-1.)/2.
  f1 = .5
  t2i = .5/edge

  for k = 0,nz-1 do begin
    zdiff1 = zc(k) - zmx1
    zdiff2 = zc(k) - zmax 
    zdiff3 = zc(k) + zmx21
    zdiff4 = zc(k) + 0.0  
    
    if(zdiff1 lt 0) then a1(k) = -zdiff1+2.*edge*alog(1.+exp( zdiff1/edge))
    if(zdiff1 ge 0) then a1(k) =  zdiff1+2.*edge*alog(1.+exp(-zdiff1/edge))
    if(zdiff2 lt 0) then a2(k) = -zdiff2+2.*edge*alog(1.+exp( zdiff2/edge))
    if(zdiff2 ge 0) then a2(k) =  zdiff2+2.*edge*alog(1.+exp(-zdiff2/edge))
    if(zdiff3 lt 0) then a3(k) = -zdiff3+2.*edge*alog(1.+exp( zdiff3/edge))
    if(zdiff3 ge 0) then a3(k) =  zdiff3+2.*edge*alog(1.+exp(-zdiff3/edge))
    if(zdiff4 lt 0) then a4(k) = -zdiff4+2.*edge*alog(1.+exp( zdiff4/edge))
    if(zdiff4 ge 0) then a4(k) =  zdiff4+2.*edge*alog(1.+exp(-zdiff4/edge))
           
    da1(k) = tanh(zdiff1*t2i) 
    da2(k) = tanh(zdiff2*t2i) 
    da3(k) = tanh(zdiff3*t2i) 
    da4(k) = tanh(zdiff4*t2i) 
     
  endfor

  n1 = tzr*tzr    + f1*(da1-da2+da3-da4)      
  t1 = tzr*tzr*zc + f1*( a1- a2+ a3- a4)      
  t1 = t1 +1.45

  e = tzr*tzr + exp(-z*z*z*z)

;# Beta them
  beta1 = (t1(nz-1) - t1(0))/(z(nz-1)-z(0))
  beta2 = (t2(nz-1) - t2(0))/(z(nz-1)-z(0))
  bint1 = t1(nz-1)-beta1*z(nz-1)
  bint2 = t2(nz-1)-beta2*z(nz-1)
  h1 = t1 - (beta1*z+bint1)
  h2 = t2 - (beta2*z+bint2)
;# plotting
;# INSERT JPEGGY STUFF
loadct,42,file='./bwah.tbl'
xpic = 1000
ypic = 800
mybgnd = fltarr(xpic,ypic)
mybgnd(*,*) = 255
mybgnd(0,0) = 0.
window,xsize=1000,ysize=600
tvscl,mybgnd(*,*),0,0

;# plot,v,z,xrange=[-.2,.8],yrange=[-27./2.,27./2.],ystyle=1,thick=2,xthick=2,ythick=2, charthick=2,charsize=2,/noerase,color=0

plot,n1,z,yrange=[-1.5,1.5],ystyle=1,thick=2,xthick=2,ythick=2,charsize=2,charthick=2,color=0,/noerase, position=[60,140,470,560],/device
oplot,n2,z,thick=2,line=2,color=0
oplot,e,z,thick=2,line=3,color=0

plot,t1,z,yrange=[-1.5,1.5],ystyle=1,thick=2,xthick=2,ythick=2,charsize=2,charthick=2,color=0,/noerase, position=[560,140,970,560],/device
oplot,t2,z,thick=2,line=2,color=0

myimage = tvrd()
help,myimage
TVLCT,r,g,b,/get
image24=bytarr(3,xpic,ypic)
image24[0,*,*]=r(myimage(*,*))
image24[1,*,*]=g(myimage(*,*))
image24[2,*,*]=b(myimage(*,*))
picfile = './new1.jpeg'
write_jpeg,picfile,image24,true=1

end
