!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! THIS IS DragoRT1.0, THE IMPLEMENTATION OF THE PBMC ALGORITHM DESCRIBED IN:
!
! PRE-CONDITIONED BACKWARD MONTE CARLO SOLUTIONS TO RADIATIVE 
! TRANSPORT IN PLANETARY ATMOSPHERES. FUNDAMENTALS: 
! SAMPLING OF PROPAGATION DIRECTIONS IN POLARISING MEDIA
! BY A. GARCIA MUÑOZ & F.P. MILLS
! ASTRONOMY & ASTROPHYSICS, XXXXXXXX, 2014
!
! THE MANUSCRIPT IS ALSO AVAILABLE AT ASTRO.PH, ID: XXXXXXXXX
!
! DragoRT1.0 M CALCULATES THE STOKES VECTOR FOR OUTGOING RADIATION 
! REFLECTED FROM A PLANE-PARALLEL ATMOSPHERE
!
! THE REFERENCE INCLUDES NUMEROUS TEST CASES THAT EXPLORE THE
! ALGORITHM'S PERFORMANCE FOR BOTH RAYLEIGH AND MIE SCATTERING MEDIA
!
! LAST MODIFIED ON 01/AUGUST/2014
!
! FOR A SKETCH OF THE RELEVANT GEOMETRICAL ANGLES, SEE FIGURE 5 IN THE REFERENCE.
! IN BRIEF: 
!    COSPOLOBS & COSPOLSUN ARE THE COSINES OF THE OBSERVER AND SOLAR POLAR ANGLES
!    AZIMUTH IS THE AZIMUTH ANGLE BETWEEN THE OBSERVER AND SOLAR PLANES.       
!
! THIS PACKAGE IS ACCOMPANIED BY:
!   * INPUT FILES: INPUT_hazeL.txt, INPUT_L13.txt, INPUT_L60.txt;
!     THEY CONTAIN EXPLANATIONS TO THE INPUT PARAMETERS. 
!     COPY INPUT_XXXX.txt INTO INPUT.dat TO EXECUTE SOME OF THE EXAMPLES DESCRIBED
!     IN THE REFERENCE
!   * FILES WITH SCATTERING MATRIX PROPERTIES: phF_hazeL.txt, phF_L13.txt, phF_L60.txt
!   * SCRIPT FOR COMPILATION IN GFORTRAN (myscript)
!
! FOR RAYLEIGH ATMOSPHERES, ONE CAN SET UP A SIMILAR phF FILE. OTHERWISE, ONE CAN
! INCORPORATE THE RELEVANT SCATTERING MATRIX BY UNCOMMENTING THE 9 ROWS FOLLOWING:
!   "FOR RAYLEIGH SCATTERING ATMOSPHERES, UNCOMMENT THE FOLLOWING 9 ROWS"
! IN THE FORTRAN CODE
!
! THE PROGRAM WILL OUTPUT THE FILE OUT.DAT WITH A SUMMARY OF THE INPUT PARAMETERS AND
! THE OUTPUT STOKES VECTORS. NOTE THAT THE DEFINITION HERE OF POSITIVE V STOKES ELEMENT 
! DIFFERS FROM THE DEFINITION IN THE SOURCES FOR THE TEST CASES.
!
! THE PROGRAM IS FREELY AVAILABLE FOR ANY NON-FOR-PROFIT APPLICATION.
! IF YOU USE IT, PLEASE GIVE CREDIT TO THE ABOVE REFERENCE.
!
! FOR COMMENTS, QUESTIONS, ETC, PLEASE CONTACT A. GARCIA MUÑOZ ON: TONHINGM@GMAIL.COM
! ALSO, IF YOU SPOT A BUG OR INCONSISTENCY, PLEASE DROP ME AN EMAIL.
! THE INTENT IS TO UPDATE DragoRT TO IMPROVE ITS PERFORMANCE AND INCLUDE NEW FEATURES,
! SO FEEL WELCOME TO DROP ME AN EMAIL TO FIND OUT ABOUT THE LATEST VERSION.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  PROGRAM PBMC
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000

  REAL, PARAMETER:: PI=4.*atan(1.)
  REAL, PARAMETER:: MACHPR=1.D-8

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  REAL IVECT0(1:4,1:NSUNMAX), HH0(1:4,1:4), weight0
  REAL IVECT0IPH(1:4,1:NSUNMAX)

  REAL X(1:3), XA(1:3), XB(1:3), SS(1:3), SP(1:3), PP(1:4,1:4)
  REAL IVECTK(1:4,1:NSUNMAX), HH(1:4,1:4), weight
  REAL IVECTCHECK

  INTEGER ISUN, JJ, IPH, IMOVE

  REAL AXXB, GEOBIAS, LL
  CHARACTER GOFOR

  INTEGER IDUM
  REAL MYRAN, EPSA, RHO
  
  CHARACTER(LEN=100) FILENAME

!!!!!!!!!!!!!

  IDUM=0
  CALL SEED(IDUM)

!!!!!!!!!!!!!
!!

  CALL SETUP

!!
!!!!!!!!!!!!!

  FILENAME='./OUT.dat' 
  OPEN(95,FILE=TRIM(FILENAME), ACCESS='APPEND', STATUS='UNKNOWN')
  WRITE(95,'(A15,20A15)') 'IPH', 'TAU', 'SSA', 'ALBGR', 'COSPOL0BS', 'COSPOLSUN', &
                          'AZIMUTH', 'STOKES I', 'STOKES Q', 'STOKES U', 'STOKES V'
  CLOSE(95)

!!!!!!!!!!!!!
!!  INITIALIZATION

  IVECT0(1:4,1:NSUNMAX)=0.

  HH0(1:4,1:4)=0.
  DO JJ=1,4
    HH0(JJ,JJ)=1.
  ENDDO

  weight0=1.

!!
!!!!!!!!!!!!!

!!!!!!!!!!!!!
!!  IPH LOOP

  IPH=0
  DO WHILE (IPH .LT. NPH)

    X=X0            
    SS=SS0 
    HH=HH0 
    weight=weight0 

!!!!!!!!!!!!!
!!  IMOVE LOOP

    IMOVE=0
    IVECT0IPH=0.
    DO WHILE (weight .GE. weightmin)

      CALL FINDXB( X, SS, XB, AXXB )

      IF (XB(1) .NE. 0. .OR. ALBGR .EQ. 0.) THEN 
        GEOBIAS=AXXB
	GOFOR  = 'A'
      ELSE
        GEOBIAS=1. 
	    RHO=MYRAN(IDUM)
          IF (AXXB .GE. RHO .AND. RHO .GT. 0.) GOFOR='A'
          IF (AXXB .LT. RHO .AND. RHO .LT. 1.) GOFOR='B'
      ENDIF

      !!!!

      IF ( GOFOR .EQ. 'A' ) THEN
        EPSA=MYRAN(IDUM)
    	LL= -1. / GAMMA * ALOG( 1. + AXXB * (-1. + EPSA) ) 
  	    XA= X -  ss * LL 
	      IF (XA(1) .LT. MACHPR)               XA(1)= 0. 
              IF (XA(1) .GT. DXSLAB*(1. - MACHPR)) XA(1)= DXSLAB

        CALL LLMMLLA( IPH, IDUM, HH, SS, SP, PP ) 
        CALL CONTRIBSUNA( XA, SS, GEOBIAS, weight, HH, IVECTk ) 
 
      !!!
 
      HH=MATMUL(HH,PP)
        HH=HH/HH(1,1)
      weight=weight * GEOBIAS * SSA
      X=XA
      SS=SP

      ELSEIF  ( GOFOR .EQ. 'B' ) THEN
        LL= -1. / GAMMA * ALOG( 1. - AXXB )
        XB= X -  ss * LL 
          IF (XB(1) .LT. MACHPR)               XB(1)= 0. 

        CALL LLMMLLB( IDUM, SP, PP ) 
        CALL CONTRIBSUNB( XB, weight, HH, IVECTk ) 

      !!!

      HH=MATMUL(HH,PP)
        HH=HH/HH(1,1)

      weight=weight * ALBGR 
      X=XB
      SS=SP

      ENDIF

      IVECT0IPH= IVECT0IPH + IVECTk

      IMOVE=IMOVE+1
    END DO

!!  IMOVE LOOP
!!!!!!!!!!!!!

    IVECTCHECK=SUM(IVECT0IPH)
    IF(ISNAN(IVECTCHECK)) THEN
      IPH=IPH
      write(87,*) IPH, 'ISNAN(IVECTkCHECK)'
    ELSE      
      IPH=IPH+1     
      IVECT0=IVECT0+IVECT0IPH    
    ENDIF

    IF ( MOD(ALOG10(FLOAT(IPH)),1.) .EQ. 0. .OR. IPH .EQ. NPH) THEN
      OPEN(95,FILE=TRIM(FILENAME), ACCESS='APPEND', STATUS='UNKNOWN')
      DO ISUN=1,NSUN
        WRITE(95,'(I15,20ES15.6)') IPH, TAU, SSA, ALBGR, COSPOL0, COSPOLSUN(ISUN), &
                                   AZIMUTH(ISUN), IVECT0(1:4,ISUN) / FLOAT(IPH)
      ENDDO
      CLOSE(95)
    ENDIF

  ENDDO

!!  IPH LOOP
!!!!!!!!!!!!!

  END PROGRAM PBMC



  SUBROUTINE SETUP
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000
  REAL, PARAMETER:: PI=4.*atan(1.)

  INTEGER NSUN, NTHETAD, NTHETAI

  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL COSPOL0, AZI0

  INTEGER ISUN

  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), AZIpSUN
  REAL SS0(1:3), SSSUN(1:3,NSUNMAX)

  INTEGER NPH
  
  REAL WEIGHTMIN, DXSLAB, TAU, SSA, ALBGR
  REAL GAMMA
  
  REAL x0(3)
  
  REAL FLUXSUN(4)
  
  INTEGER ITHETAD, ITHETAI
  
  REAL THETAD(0:NTHETADMAX), DTHETAD, MUD(0:NTHETADMAX)
  
  REAL A1(0:NTHETADMAX), B1(0:NTHETADMAX), A3(0:NTHETADMAX), B2(0:NTHETADMAX) 
  REAL MM0(4,4,0:NTHETADMAX)
  
  REAL FFD(0:NTHETADMAX),FFI(0:NTHETAIMAX),DFFI,THETAI(0:NTHETAIMAX)
  
  CHARACTER(LEN=20) FILESCATMATRIX

  DATA ex / 1.,0.,0. /, ey / 0.,1.,0. /, ez / 0.,0.,1. /
  DATA FLUXSUN / pi,0.,0.,0./

  OPEN(99,FILE='./INPUT.dat')
    READ(99,*) COSPOL0
      AZI0=00.              ! TO HAVE OBSERVER IN THE XZ PLANE
      IF(COSPOL0 .EQ. 1.) COSPOL0=0.99999999D0

    SS0= COSPOL0 * ex +                                                             &
         SQRT(1.-COSPOL0**2.) * ( COS(AZI0*pi/180D0) * ez - SIN(AZI0*pi/180D0) * ey ) 

    READ(99,*) NSUN
!      IF (NSUN .NE. NSUNMAX) THEN
!        WRITE(*,*) 'NSUN .NE. NSUNMAX'
!	    STOP
!      ENDIF
    DO ISUN=1,NSUN
      READ(99,*) COSPOLSUN(ISUN), AZIMUTH(ISUN)
      AZIpSUN= 180D0 + AZI0 - AZIMUTH(ISUN)
      SSSUN(1:3,ISUN)= COSPOLSUN(ISUN) *ex +                                          &
      SQRT(1.-COSPOLSUN(ISUN)**2.) * ( COS(AZIpSUN*pi/180D0)*ez - SIN(AZIpSUN*pi/180D0)*ey )
    ENDDO
    SSSUN=-SSSUN

    READ(99,*) NPH
    READ(99,*) WEIGHTMIN
    READ(99,*) DXSLAB
    READ(99,*) TAU
    READ(99,*) SSA
    READ(99,*) ALBGR
    
    READ(99,*) FILESCATMATRIX

    x0(1)  = DXSLAB
    x0(2:3)= 0.
    
    GAMMA=TAU / DXSLAB
    
    READ(99,*) NTHETAD
      IF (NTHETAD .NE. NTHETADMAX) THEN
        WRITE(*,*) 'NTHETAD .NE. NTHETADMAX'
	    STOP
      ENDIF

     DTHETAD= pi / FLOAT(NTHETAD)
     OPEN(98,FILE=FILESCATMATRIX)
     DO ITHETAD=NTHETAD,0,-1
       READ(98,*) THETAD(ITHETAD),&
       A1(ITHETAD),A3(ITHETAD),B1(ITHETAD),B2(ITHETAD)
     ENDDO
     CLOSE(98)
     THETAD=THETAD * pi /180.
     MUD=COS(THETAD)

!    FOR RAYLEIGH SCATTERING ATMOSPHERES, UNCOMMENT THE FOLLOWING 9 ROWS
!    DTHETAD= pi / FLOAT(NTHETAD)
!    DO ITHETAD=0,NTHETAD
!      THETAD(ITHETAD)=pi - FLOAT(ITHETAD) * DTHETAD
!    ENDDO 
!    MUD=COS(THETAD)
!    A1(:)=3./4. * ( 1. + MUD(:)**2.)
!    B1(:)=3./4. * (-1. + MUD(:)**2.) 
!    A3(:)=3./2. * MUD(:) 
!    B2(:)=0.  

    MM0=0.
    MM0(1,1,:)= A1
    MM0(1,2,:)= B1
    MM0(2,1,:)= B1
    MM0(2,2,:)= A1
    MM0(3,3,:)= A3 
    MM0(3,4,:)= B2
    MM0(4,3,:)=-B2  
    MM0(4,4,:)= A3
      MM0=MM0 / 4. / PI

    READ(99,*) NTHETAI
      IF (NTHETAI .NE. NTHETAIMAX) THEN
        WRITE(*,*) 'NTHETAI .NE. NTHETAIMAX'
	    STOP
      ENDIF

  CLOSE(99)

  CALL FFDGENERATOR( A1, NTHETAD, THETAD, DTHETAD, FFD) ! PRODUCES FFD
  CALL THETAIGENERATOR( NTHETAD, THETAD, FFD, NTHETAI, FFI, DFFI, THETAI ) ! PRODUCES THETAI

  END SUBROUTINE SETUP


  SUBROUTINE FFDGENERATOR( A1, NTHETAD, THETAD, DTHETAD, FFD)  ! EXTENDED TRAPEZOIDAL RULE
  implicit none

  INTEGER ITHETAD, NTHETAD
  REAL A1(0:NTHETAD), THETAD(0:NTHETAD), FFD(0:NTHETAD), DTHETAD

  FFD(0)=0.
  DO ITHETAD=1, NTHETAD
  FFD(ITHETAD)= FFD(ITHETAD-1) + &
    ( A1(ITHETAD-1) * SIN(THETAD(ITHETAD-1)) + A1(ITHETAD) * SIN(THETAD(ITHETAD)) ) 
  ENDDO  

  FFD= FFD * DTHETAD / 2. / 2.
  FFD(NTHETAD)=1.

  END SUBROUTINE FFDGENERATOR


  SUBROUTINE THETAIGENERATOR( NTHETAD, THETAD, FFD, NTHETAI, FFI, DFFI, THETAI )
  implicit none

  INTEGER NTHETAD, NTHETAI, ITHETAI
  REAL FFD(0:NTHETAD), THETAD(0:NTHETAD), DER2(0:NTHETAD), FFI(0:NTHETAI), DFFI, THETAI(0:NTHETAI), DUMMY

    DFFI=1./FLOAT(NTHETAI)    
    DO ITHETAI=0,NTHETAI
      FFI(ITHETAI)= FLOAT(ITHETAI) * DFFI 
    ENDDO

    CALL SPLINE(FFD(0:NTHETAD),THETAD(0:NTHETAD),NTHETAD+1,2.E30,2.E30,DER2(0:NTHETAD))
 
    DO ITHETAI=0,NTHETAI
      CALL SPLINT(FFD(0:NTHETAD),THETAD(0:NTHETAD),DER2(0:NTHETAD),NTHETAD+1, &
                  FFI(ITHETAI),THETAI(ITHETAI),DUMMY)
    ENDDO

  END SUBROUTINE THETAIGENERATOR


  SUBROUTINE FINDXB( X, SS, XB, AXXB )
  implicit none
  
  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI


  REAL X(1:3), SS(1:3), XB(1:3), AXXB
  REAL MUX, LL

    MUX=SS(1)

    IF (MUX .GT. 0.) THEN       ! -SS intersects the surface
      LL=(X(1)-    0.)/MUX 
      XB= X-SS*LL 
      XB(1)=0.
    ELSEIF (MUX .LT. 0.) THEN   ! -SS intersects the top of atmosphere
      LL=(X(1)-DXSLAB)/MUX 
      XB= X-SS*LL 
      XB(1)=DXSLAB
    ELSEIF (MUX .EQ. 0.) THEN   ! -SS is parallel to yz plane
      WRITE(*,*) 'MUX = 0 @ FINDXB' 
      STOP
    ENDIF
   
    AXXB=1. - exp(-LL*GAMMA)
      
  END SUBROUTINE FINDXB



  SUBROUTINE LLMMLLA( IPH, IDUM, HH, SS, SP, PP ) 
  implicit none
  
  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000

  INTEGER NPH, NSUN, NTHETAD, NTHETAI, IPH
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  INTEGER IDUM
  REAL HH(1:4,1:4), PP(1:4,1:4), SS(1:3), SP(1:3)
  REAL MM(1:4,1:4)
  
  REAL MYRAN
  REAL ETAA, THETAA, COSTHETAA
  
  REAL ZETAA, YZETAA, PHIA, gg
  REAL m11, m12, qq, uu
  INTEGER REJ

  REAL, PARAMETER:: PI=4.*atan(1.)
  REAL, PARAMETER:: MACHPR=1.D-8

  REAL e1(1:3), e2(1:3), e3(1:3)
  REAL ii, ip
  REAL NORM, DOTP

  REAL COSAA, COSAP, SINAP, THETACORR, COSIP
  REAL LLkk(1:4,1:4), LLkp(1:4,1:4)

  INTEGER INDL, INDR

  ETAA=MYRAN(IDUM)  
    IF(ETAA .EQ. 0.) ETAA=MYRAN(IDUM)
    IF(ETAA .EQ. 1.) ETAA=MYRAN(IDUM)
  CALL FINDTHETAA(ETAA, NTHETAI, THETAI, FFI, DFFI, THETAA)  !! INVERSE: GIVEN ETAA --> FIND THETAA=THETAI(FFI)

  CALL MMMATRIX(THETAA, MM)


  !!!

  m11=MM(1,1)
  m12=MM(1,2) 
  qq=HH(1,2)/HH(1,1) 
  uu=HH(1,3)/HH(1,1)

  REJ=1
  DO WHILE (REJ .EQ. 1)
  ZETAA=MYRAN(IDUM)
  PHIA =2. * PI * ZETAA
  
  YZETAA=2. * MYRAN(IDUM)
  gg= 1. + m12/m11 * ( qq * cos(2.*PHIA) - uu * sin(2.*PHIA) )

  IF (YZETAA .LE. gg) REJ=0
  END DO


  IF     (PHIA .LT. PI) THEN
    ii= PI - PHIA   
  ELSEIF (PHIA .GE. PI) THEN 
    ii= 2. * PI - PHIA
  ENDIF

  e3=SS/NORM(SS)

  CALL CROSSP(e3,ez,e2)
  e2=e2/NORM(e2)

  CALL CROSSP(e2,e3,e1)
  e1=e1/NORM(e1)

  !!!

  SP=COS(THETAA)*e3 + SIN(THETAA) * (COS(PHIA)*e1 + SIN(PHIA)*e2)
  SP=SP/NORM(SP)

  !!!

  COSAA=DOTP(SS,ez)  
  COSAP=DOTP(SP,ez) 
  SINAP=SQRT(1. - COSAP**2.)

  IF     (PHIA .LT. PI) THEN
    THETACORR= 2.* PI - THETAA 
  ELSEIF (PHIA .GE. PI) THEN 
    THETACORR=THETAA
  ENDIF

  COSIP=(COSAA - COSAP * COS(THETACORR)) / (SINAP * SIN(THETACORR)) 
  IF(COSIP .GE. 1.) THEN
    write(87,*) 'COSIP GE 1', COSIP
    COSIP=1.
!    pause
  ELSEIF(COSIP .LE. -1.) THEN  
    write(87,*) 'COSIP LE -1', COSIP  
    COSIP=-1.
!    pause
  ENDIF
  
  IP=ACOS(COSIP)

  CALL ROTATION (LLkk, PI-ii)
  CALL ROTATION (LLkp,   -ip)
  
  PP=MATMUL(LLkk,MM)
  PP=MATMUL(PP,LLkp) 

  END SUBROUTINE LLMMLLA


  SUBROUTINE LLMMLLB( IDUM, SP, PP ) 
  implicit none
  
  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  INTEGER IDUM
  REAL PP(1:4,1:4), SP(1:3)

  REAL MYRAN
  REAL ETAB, THETAB, COSTHETAB

  REAL ZETAB, PHIB
  REAL NORM
  
  REAL, PARAMETER:: PI=4.*atan(1.)

  ETAB=MYRAN(IDUM)  
  COSTHETAB=SQRT(ETAB)
  THETAB=ACOS(COSTHETAB)

  !!!

  ZETAB=MYRAN(IDUM)
  PHIB =2. * PI * ZETAB

  !!!

  SP=COS(THETAB)*ex + SIN(THETAB) * (COS(PHIB)*ez + SIN(PHIB)*ey)
  SP=-SP/NORM(SP)

  !!!
  
  PP=0.
  PP(1,1)=1.

  END SUBROUTINE LLMMLLB



  SUBROUTINE FINDTHETAA(ETAA, NTHETAI, THETAI, FFI, DFFI, THETAA)
  implicit none
  
  INTEGER NTHETAI, INDL, INDR
  REAL ETAA, THETAI(0:NTHETAI), FFI(0:NTHETAI), DFFI, THETAA, COSTHETAA
  
  REAL, PARAMETER:: MACHPR=1.D-8
  REAL, PARAMETER:: PI=4.*atan(1.)
  
  INDL=INT(ETAA/DFFI)
  INDR=INDL+1
  COSTHETAA=COS(THETAI(INDL)) + &
        (COS(THETAI(INDR))-COS(THETAI(INDL))) / (FFI(INDR)-FFI(INDL)) * &
	(ETAA-FFI(INDL))
  THETAA=ACOS(COSTHETAA)	

  IF(THETAA .LE. 0.) THETAA=PI* MACHPR
  IF(THETAA .GE. PI) THETAA=PI*(1.-MACHPR)

  END SUBROUTINE FINDTHETAA


  SUBROUTINE MMMATRIX(THETA, MM)
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000
  
  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  REAL, PARAMETER:: PI=4.*atan(1.)

  INTEGER INDL, INDR
  REAL MM(1:4,1:4), THETA

  INDL=NTHETAD-INT(THETA/DTHETAD)
  INDR=INDL-1

  MM(1:4,1:4)= MM0(1:4,1:4,INDL) + &
    (MM0(1:4,1:4,INDR)-MM0(1:4,1:4,INDL))/(COS(THETAD(INDR))-COS(THETAD(INDL))) * &
    (COS(THETA)-COS(THETAD(INDL)))

  END SUBROUTINE MMMATRIX



  SUBROUTINE CONTRIBSUNA( X, SS, GEOBIAS, weight, HH, IVECTk )
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000
  
  REAL, PARAMETER:: PI=4.*atan(1.)
  REAL, PARAMETER:: MACHPR=1.D-8

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  INTEGER ISUN
  REAL X(1:3), SS(1:3), GEOBIAS, weight, HH(1:4,1:4), IVECTk(1:4,1:NSUNMAX)
  REAL XSUN(1:3), AXXSUN, TXXSUN, iiSUN, ipSUN
  REAL LLkkSUN(1:4,1:4), LLkpSUN(1:4,1:4), MMSUN(1:4,1:4), PPSUN(1:4,1:4)
  REAL DOTP

  DO ISUN=1,NSUN

  CALL FINDXB( X, SSSUN(1:3,ISUN), XSUN, AXXSUN )  
  TXXSUN=1. - AXXSUN
  
  CALL METHOD2( SS, SSSUN(1:3,ISUN), iiSUN, ipSUN )

  CALL ROTATION(LLkkSUN, PI-iiSUN)
  CALL ROTATION(LLkpSUN,   -ipSUN)

  CALL MMMATRIX(ACOS(DOTP(SS,SSSUN(1:3,ISUN))),MMSUN)

  PPSUN=MATMUL(LLkkSUN,MMSUN)
  PPSUN=MATMUL(PPSUN,LLkpSUN)

  IVECTk(1:4,ISUN)= weight * MATMUL(MATMUL(HH,PPSUN), FLUXSUN ) *  TXXSUN * SSA * GEOBIAS
  
  ENDDO

  END SUBROUTINE CONTRIBSUNA


  SUBROUTINE CONTRIBSUNB( X, weight, HH, IVECTk )
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000
  
  REAL, PARAMETER:: PI=4.*atan(1.)
  REAL, PARAMETER:: MACHPR=1.D-8

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI

  INTEGER ISUN
  REAL X(1:3), weight, HH(1:4,1:4), IVECTk(1:4,1:NSUNMAX)
  REAL XSUN(1:3), AXXSUN, TXXSUN
  REAL DOTP

  DO ISUN=1,NSUN

  CALL FINDXB( X, SSSUN(1:3,ISUN), XSUN, AXXSUN )  
  TXXSUN=1. - AXXSUN

  IVECTk(1:4,ISUN)= weight * MATMUL(HH,FLUXSUN) *  TXXSUN * ALBGR / PI * DOTP(-ex,SSSUN(1:3,ISUN))
  
  ENDDO
  
  END SUBROUTINE CONTRIBSUNB


  SUBROUTINE METHOD2( SS, SP, ii, ip )
  implicit none

  COMMON /PARAMETERS1/ NPH, NSUN, NTHETAD, NTHETAI, &
         ex, ey, ez, WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB, &
         X0, COSPOL0, SS0, COSPOLSUN, AZIMUTH, SSSUN, FLUXSUN, &
	     MM0, THETAD, DTHETAD, THETAI, FFI, DFFI

  INTEGER, PARAMETER:: NSUNMAX=10, NTHETADMAX= 1800, NTHETAIMAX=4000

  INTEGER NPH, NSUN, NTHETAD, NTHETAI
  REAL ex(1:3), ey(1:3), ez(1:3)
  REAL WEIGHTMIN, TAU, GAMMA, SSA, ALBGR, DXSLAB
  REAL X0(1:3), COSPOL0, SS0(1:3)
  REAL COSPOLSUN(1:NSUNMAX), AZIMUTH(1:NSUNMAX), SSSUN(1:3,1:NSUNMAX), FLUXSUN(1:4)
  REAL MM0(1:4,1:4,0:NTHETADMAX), THETAD(0:NTHETADMAX), DTHETAD
  REAL THETAI(0:NTHETAIMAX), FFI(0:NTHETAIMAX), DFFI
    
  REAL SS(1:3), SP(1:3), ii, ip
  REAL VMERID(1:3), VSCATT(1:3)
  
  REAL DOTP, NORM  
  

  CALL CROSSP(SS,ez,VMERID)
  VMERID=VMERID/NORM(VMERID)
  
  CALL CROSSP(SS,SP,VSCATT)
  VSCATT=VSCATT/NORM(VSCATT)
  
  IF (DOTP(VSCATT,ez) .LT. 0.) VSCATT=-VSCATT
  ii=ACOS(DOTP(VMERID,VSCATT))  

  !!!

  CALL CROSSP(SP,ez,VMERID)
  VMERID=VMERID/NORM(VMERID)

  CALL CROSSP(SP,SS,VSCATT)
  VSCATT=VSCATT/NORM(VSCATT)

  IF (DOTP(VSCATT,ez) .GT. 0.) VSCATT=-VSCATT
  ip=ACOS(DOTP(VMERID,VSCATT))
  
  END SUBROUTINE METHOD2


  SUBROUTINE seed(idum)
  implicit none
  
  integer idum ,i,da_ti_zo(8)
  character(8) date
  character(10) time
  character(5) zone
  
  idum=0
  
  call date_and_time(date,time,zone,da_ti_zo)

  do i=1,8
   idum=idum+da_ti_zo(i)
  enddo
  
  idum=-idum  

!  idum=-1 ! THIS LINE MUST BE COMMENTED TO ALLOW THE SEEED TO EFFECTIVELY CHANGE OVER TIME
  open(99,file='seed.out') 
    write(99,*) 'SEED EFFECTIVELY USED IN THE CALCULATIONS=', idum
    write(*,*)  ''
    write(*,*)  'SEED EFFECTIVELY USED IN THE CALCULATIONS=', idum
!    write(*,*)  'LINE ABOVE WITH idum=-1 SHOULD EVENTUALLY BE COMMENTED'
    write(*,*)  ''    
  close(99)
  

  END SUBROUTINE seed


  SUBROUTINE CROSSP (A, B, C)                                            
  implicit none

  REAL A(1:3), B(1:3), C(1:3)

  C(1) = A(2)*B(3) - A(3)*B(2)                                              
  C(2) = A(3)*B(1) - A(1)*B(3)
  C(3) = A(1)*B(2) - A(2)*B(1)

  RETURN
  END SUBROUTINE CROSSP


  SUBROUTINE ROTATION (LL, ii)                                            
  implicit none

  REAL LL(1:4,1:4), ii

    LL(1:4,1:4)=0.
  
    LL(1,1)=1.
  
    LL(2,2)=cos(2.*ii)
    LL(3,3)=LL(2,2)
    LL(2,3)=sin(2.*ii)
    LL(3,2)=-LL(2,3)
  
    LL(4,4)=1.

  RETURN
  END SUBROUTINE ROTATION

      
  REAL FUNCTION NORM(a)
  implicit none
  
  REAL a(1:3)
  
    NORM= SQRT( a(1)**2 + a(2)**2 + a(3)**2 ) 
  
  RETURN
  END FUNCTION NORM


  REAL FUNCTION DOTP(a,b)
  implicit none

  REAL a(1:3), b(1:3)

    DOTP= a(1)*b(1) + a(2)*b(2) + a(3)*b(3) 
  
  RETURN
  END FUNCTION DOTP


  SUBROUTINE spline(x,y,n,yp1,ypn,y2) ! FROM NUMERICAL RECIPES IN FORTRAN 77. PRESS ET AL. ISBN 0 521 43064 X
  INTEGER n,NMAX
  REAL yp1,ypn,x(n),y(n),y2(n)
  PARAMETER (NMAX=500000)
  !! Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y_i=f(x_i)m with
  !! x1<x2<...<xN, and given values yp1 and ypn for the first derivative of the inter-
  !! polating function at points 1 and n, respectively, this routine returns an array y2(1:n) of
  !! length n which contains the second derivatives of the interpolating function at the tabulated
  !! points x_i. If yp1 and/or ypn are equal to 1E30 or larger, the routine is signaled to set
  !! the corresponding boundary condition for a natural spline, with zero second derivative on
  !! that boundary
  !! Parameter NMAX is the largest anticipated value of n.
  INTEGER i,k
  REAL p,qn,sig,un,u(NMAX)
  if (yp1.gt..99e30) then  !the lower boundary condition is set either to be "natural"
      y2(1)=0.  
      u(1)= 0.
  else
      y2(1)=-0.5
      u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
  endif
  do i=2,n-1    !this is the decomposition loop of the tridiagonal algorithm. y2 and u are used
    sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))         !for temporary storage of the decomposed factors
    p=sig*y2(i-1)+2.
    y2(i)=(sig-1.)/p
    u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))  &
&        /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
  enddo
  if (ypn.gt..99e30) then  !the upper boundary condition is set either to be "natural"
      qn=0.
      un=0.
  else
      qn=0.5
      un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
  endif
  y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
  do k=n-1,1,-1
       y2(k)=y2(k)*y2(k+1)+u(k)
  enddo 
!  return
  END


  SUBROUTINE splint(xa,ya,y2a,n,x,y,dydx)  ! FROM NUMERICAL RECIPES IN FORTRAN 77. PRESS ET AL. ISBN 0 521 43064 X
  INTEGER n
  REAL x,y,dydx,xa(n),y2a(n),ya(n)
!! Given te arrays xa(1:n) and ya(1:n) of length n, which tabulate a function (with the
!! xa_i's in order), and given the array y2a(1:n), which is the output from spline above,
!! and given a value of x this routine returns a cubic-spline interpolated value y.
  INTEGER k,khi,klo
  REAL a,b,h
  klo=1       !We will find the right place in the table by means of bisection.
  khi=n       !This is optimal if sequential calls to this routine are at random
1 if(khi-klo.gt.1) then
  k=(khi+klo)/2
  if(xa(k).gt.x) then
     khi=k
  else
     klo=k
  endif
  goto 1
  endif
  h=xa(khi)-xa(klo)
  if (h.eq.0.) pause 'bad xa input in splint'  !The xa's must be distinct.
  a=(xa(khi)-x)/h
  b=(x-xa(klo))/h
  y=a*ya(klo)+b*ya(khi)+  & 
&         ((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**2)/6.   
!  dydx=(-ya(klo)+ya(khi)+(h**2)/6.*(-(3*aa**2-1)*y2a(klo)+(3*bb**2-1)*y2a(khi)))/h
  return
  END
      
  FUNCTION MYRAN(idum)  ! FROM NUMERICAL RECIPES IN FORTRAN 77. PRESS ET AL. ISBN 0 521 43064 X
  !Returns a uniform random deviate between 0.0 and 1.0. Set idum to any negative value
  !to initialize or reinitialize the sequence.
  INTEGER idum
  INTEGER MBIG,MSEED,MZ
!C REAL MBIG,MSEED,MZ
  REAL myran,FAC
  PARAMETER (MBIG=1000000000,MSEED=161803398,MZ=0,FAC=1./MBIG)
!C PARAMETER (MBIG=4000000.,MSEED=1618033.,MZ=0.,FAC=1./MBIG)
  !According to Knuth, any large mbig, and any smaller (but still large) mseed can be substituted
  !for the above values.
  INTEGER i,iff,ii,inext,inextp,k
  INTEGER mj,mk,ma(55) !The value 55 is special and should not be modified; see
!C REAL mj,mk,ma(55) !Knuth.
  SAVE iff,inext,inextp,ma
  DATA iff /0/
  if(idum.lt.0.or.iff.eq.0)then !Initialization.
  iff=1
  mj=abs(MSEED-abs(idum)) !Initialize ma(55) using the seed idum and the large nummj=
  mj=mod(mj,MBIG)            !ber mseed.
  ma(55)=mj
  mk=1
  do i=1,54 !Now initialize the rest of the table,
    ii=mod(21*i,55) !in a slightly random order,
    ma(ii)=mk !with numbers that are not especially random.
    mk=mj-mk
    if(mk.lt.MZ)mk=mk+MBIG
                !274 Chapter 7. Random Numbers
                !Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X)
                !Copyright (C) 1986-1992 by Cambridge University Press. Programs Copyright (C) 1986-1992 by Numerical Recipes Software.
                !Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copying of machinereadable
                !files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website
                !http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America).
    mj=ma(ii)
  enddo
  do k=1,4 !We randomize them by \u201cwarming up the generator.\u201d
  do i=1,55
    ma(i)=ma(i)-ma(1+mod(i+30,55))
    if(ma(i).lt.MZ)ma(i)=ma(i)+MBIG
  enddo
  enddo
  inext=0 !Prepare indices for our first generated number.
  inextp=31 !The constant 31 is special; see Knuth.
  idum=1
  endif
  inext=inext+1 !Here is where we start, except on initialization. Increment
  if(inext.eq.56)inext=1 !inext, wrapping around 56 to 1.
  inextp=inextp+1 !Ditto for inextp.
  if(inextp.eq.56)inextp=1
  mj=ma(inext)-ma(inextp) !Now generate a new random number subtractively.
  if(mj.lt.MZ)mj=mj+MBIG !Be sure that it is in range.
  ma(inext)=mj !Store it,
  myran=mj*FAC !and output the derived uniform deviate.
  return
  END
