program isochron
use sep
implicit none

real			:: isoch(601),x(601),dzdx(601),x0(601),l(601)
real			:: a,b,a2,b2,vel,h,th
integer			:: i

call initpar()
call from_param("h",h,1.)
call from_param("th",th,1.)

vel=2.
a=(vel*th/2)
a2=a*a
b2=a2-h*h
b=sqrt(b2)
if(b==0) call seperr('For the specified values of h,th,vel the isochron is zero!!!')

isoch=0
do i=-300,300
   x(i+301)=i*0.01
   if (a2 >= x(i+301)*x(i+301)) then
      isoch(i+301)=-b/a*sqrt(a2-x(i+301)*x(i+301))
   end if
end do

dzdx=0
where(isoch .ne.0 )
   dzdx=-b/a*x/abs(isoch)
end where

x0=0.
where(dzdx .ne.0 )
   x0=abs(isoch)*dzdx+x
end where

l=0
where(x0.ne.0 )
   l=2*sqrt((x-x0)**2+isoch**2)
end where

l(301)=2*abs(isoch(301))

call to_history("n1",601); call to_history("d1",0.01); call to_history("o1",-3)  
call to_history("n2",1); call to_history("d2",1.); call to_history("o2",0.)
call to_history("n3",1); call to_history("d3",1.); call to_history("o3",0.)
call to_history("esize",4)
call srite("out",isoch,4*601)

call to_history("n1",601,"t0"); call to_history("d1",0.01,"t0"); call to_history("o1",-3,"t0")  
call to_history("n2",1); call to_history("d2",1.,"t0"); call to_history("o2",0.,"t0")
call to_history("n3",1,"t0"); call to_history("d3",1.,"t0"); call to_history("o3",0.,"t0")
call to_history("esize",4)
call srite("t0",l,4*601)

end program
