      PROGRAM ue5b
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION R(1000),HIST(1001)
      CHARACTER*10 TX
      CHARACTER*5 TY
      CHARACTER*40 CAPT
C identify program to user
      WRITE(*,*)' Program ue5b simulates (1+alpha*cos^2theta) distr.'
      WRITE(*,*)' Inversion Methode' 
      WRITE(*,*)' '
C ask for numerical input
      WRITE(*,*)' Enter number NEXP of trials'
      WRITE(*,*)' >'
      READ(*,*)NEXP
      WRITE(*,*)' Enter parameter alpha'
      WRITE(*,*)' >'
      READ(*,*) alpha
C setup distribution function
      beta= 0.5/(1.+alpha/3.)
      gamma= 1.+alpha/3.
      n=1
C initialization
      CAPT='alpha = ******, N_exp# = *******'
      WRITE(CAPT(9:14),'(F6.3)')alpha
      WRITE(CAPT(25:32),'(I7)')NEXP
      X0=-1.0D0
      DELX=0.1D0
      NX=20
      CALL SMHSIN(HIST,X0,DELX,NX)
      pi=atan2(0.,-1.)
c      write(*,*)pi
C simulate trials and fill histogram
      DO 20 I=1,NEXP
        CALL RNECUY(R,N)
        q= ( 3.*(gamma-r(1)/beta)/alpha/2.)
        p=1./alpha
        d=p**3 + q**2
c unterscheide d>0, d<0
        if( d.lt.0.) then
           sqd=sqrt(-d)
           a=sqrt(-d+q**2)
           aa=2.*a**(1./3.)
           phi=atan2(sqd,-q)
           cost=aa*cos((phi+pi)/3.)
c die beiden anderen Loesungen sind >+1: aa*cos(phi/3.); aa*cos((phi-pi)/3.) 
 26        continue
        else
           sqd=sqrt(d)
           s1=sign(1.,-q+sqd)*(abs(-q+sqd))**(1./3.)
           s2=sign(1.,-q-sqd)*(abs(-q-sqd))**(1./3.)
           cost=s1+s2 
        endif
        CALL SMHSFL(HIST,X0,DELX,NX,DBLE(cost),1.D0)
20    CONTINUE
C graphical output
      CALL GRNBWS()
      WRITE(*,*)' Please, enter number of workstation:'
      WRITE(*,*)' >'
      READ(*,*)NWS
C present histogram
      TX='cos(theta)'
      TY='N'
      CALL SMHSGR(HIST,X0,DELX,NX,TX,10,TY,5,CAPT,40,NWS)
      END
