program find_alpha

print*,'Enter phi and phib'
read*,phi,phib

phi=phi/45.*atan(1.)
phib=phib/45.*atan(1.)

psib=asin(sin(phib)/sin(phi))/2.-phib/2.

sig=1.
nint=100
start=-phib
end=phib
n=0

1 n=n+1
do i=1,nint
alpha=start+(end-start)*float(i-1)/(nint-1)
psi0=asin(sin(alpha)/sin(phi))/2.-alpha/2.
beta=psib-psi0-alpha
if (beta*sig.lt.0.) then
start=alpha-(end-start)/(nint-1)
end=alpha
if (n.le.3) goto 1
print*,beta,alpha*45./atan(1.)
goto 2
endif
enddo

2 continue

sig=-1.
nint=100
start=phib
end=-phib
n=0

3 n=n+1
do i=1,nint
alpha=start+(end-start)*float(i-1)/(nint-1)
psi0=asin(sin(alpha)/sin(phi))/2.-alpha/2.
beta=psib-psi0-alpha
if (beta*sig.lt.0.) then
start=alpha-(end-start)/(nint-1)
end=alpha
if (n.le.3) goto 3
print*,beta,alpha*45./atan(1.)
stop
endif
enddo

end
