!***************************************************** !* Program to demonstrate the hyperbolic subroutines * !* ------------------------------------------------- * !* Reference: BASIC Scientific Subroutines, Vol. II * !* By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * !* * !* F90 version by J-P Moreau. * !* (www.jpmoreau.fr) * !***************************************************** Program Hyper integer ligne real*8 XX,Y real*8 CosH,SinH,TanH character*1 r ligne=1 print *,' ' print *,' X SINH(X) COSH(X) TANH(X) ' print *,' ---------------------------------------------------' XX=-5.d0 do while (XX<=5.2d0) write(*,10,advance='no') XX Y=SinH(XX) write(*,11,advance='no') Y Y=CosH(XX); write(*,11,advance='no') Y Y=TanH(XX); write(*,11,advance='no') Y if (ligne==20) then ligne=0 read *, r else print *,' ' end if ligne=ligne+1 XX=XX+0.2d0 end do print *,' ' print *,' ' stop 10 format(F6.1,' ') 11 format(F14.10,' ') end !------------------------------------------------------------ ! Modified cordic exponential subroutine ! ---------------------------------------------------------- ! This subroutine takes an input value and returns Y:=EXP(X) ! X may be any positive or negative real value !------------------------------------------------------------ real*8 Function XExp(X) COMMON /C1/ E,Z real*8 X, E,K,Y,Z, A(1:9),W(1:9) integer I ! Get coefficients call Coeffs(N,A) ! Reduce the range of X K=int(X) X=X-K ! Determine the weighting coeffs, W(I) call WCoeffs(N,X,W) ! Calculate products Y=1.d0 do I=1, N if (W(I)>0.d0) Y=Y*A(I) end do ! Perform residual multiplication Y=Y*(1.d0+Z*(1.d0+Z/2.d0*(1.d0+Z/3.d0*(1.d0+Z/4.d0)))) ! Account for factor EXP(K) if (K<0.d0) E=1.d0/E if (dabs(K)<1.d0) then XExp=Y return end if do I=1, INT(dabs(K)) Y=Y*E end do ! Restore X X=X+K XExp=Y return end Subroutine WCoeffs(N,X,W) COMMON /C1/ E,Z real*8 X, aa,E,Z, W(1:9) integer I aa=0.5d0 Z=X do I=1, N W(I)=0.d0 if (Z>aa) W(I)=1.d0 Z=Z-W(I)*aa aa=aa/2.d0 end do return end Subroutine Coeffs(N,A) COMMON /C1/ E,Z real*8 E, A(1:9) N=9 E=2.718281828459045 A(1)=1.648721270700128 A(2)=1.284025416687742 A(3)=1.133148453066826 A(4)=1.064494458917859 A(5)=1.031743407499103 A(6)=1.015747708586686 A(7)=1.007843097206448 A(8)=1.003913889338348 A(9)=1.001955033591003 return end !*---------------------------------------------* !* Hyperbolic sine Function * !* ------------------------------------------- * !* This Procedure uses the definition of the * !* hyperbolic sine and the modified cordic * !* exponential Function XExp(X) to approximate * !* SINH(X) over the entire range of real X. * !*---------------------------------------------* real*8 Function SinH(X) !Label: 10 real*8 X,Y,Z; integer I ! Is X small enough to cause round off erroe ? if (dabs(X)<0.35d0) goto 10 ! Calculate SINH(X) using exponential definition Y=XExp(X) Y=(Y-(1.d0/Y))/2.d0 SinH=Y return ! series approximation (for X small) 10 Z=1.d0; Y=1.d0 do I=1, 8 Z=Z*X*X/((2*I)*(2*I+1)) Y=Y+Z end do Y=Y*X SinH=Y return end ! hyperbolic cosine Function real*8 Function CosH(X) real*8 X,Y Y=XExp(X) Y=(Y+(1.d0/Y))/2.d0 CosH=Y return end ! hyperbolic tangent Function ! TANH(X]:=SINH(X)/COSH(X) real*8 Function TanH(X) real*8 X,V,Y V=SinH(X) Y=CosH(X) TanH=V/Y return end ! End of file hyper.f90