c-----------------------------------------------------------
c     mcnewd_sub.for
c     19.12.2001 
c-----------------------------------------------------------
c     This file contains subroutines for calculation of scattering
c     from a given particle field illuminated from a single direction,
c     mcnewd_sub.for: comments for Luke ( <--mcnew_sub_output2.for)     

c-------------------------------------------------------------
c     Main programs using this subroutine:

c      mcnewd_asyobs.for - asymmetry driver 
c                          uses no files created in this subroutine 

c      mcnewd_oppo.for   - opposition curve driver   
c                          uses wwwfile
                         
c      mcnewd_tilt.for   - tilt curve driver   
c                          uses no files created in this subroutine 

c-------------------------------------------------------------
c     Include files:
c      mcnew_sub.inc  = commons for transfer between subroutines
c      mcnew_main.inc = commons for transfer between main/subroutines
c                       includes variables related to Saturn's illumination

c-------------------------------------------------------------
c     Variables transferred via call to main program
c-------------------------------------------------------------

c     general input:

c      nphot   = number of photons
c      posfile = name of ascii-file containing particle positions,
c                particle radii, and box-dimensions
c      parfile = name of file containing some additional input params:
c         ssca    - single scattering albedo
c         maxscat - max number of scatterings
c         nbox    - max number of replica regions in search of scattering
c         nboxvis - max number of replica in visibility test

c      sunfii  = sun longitude (degrees, ccw from x-axis)
c      sunlat  = sun latitude  (degrees)

c      idum    = seed for random numbers
c      icheck  = whether or not to output some debug info
c      method  = phase-function model, explained below
c      ppar1,ppar2,ppar3 = parameters of phase function

c      mx,my   = number of divisions in spatial partitioning
c      sredu   = reduce particle sizes by this factor (SET THIS TO 1)
c      zredu   = multiply partcle z-coordinates by this factor (SET THIS TO 1)

c     input related to indirect MC
c      ntest          = number of viewing directions
c      testfii(ntest) = observing longitudes (degrees, ccw from x-axis)
c      testlat(ntest) = observing latitudes  (degrees)

c-------------------------------------------------------------
c     output related to indirect MC 
c      testw(ntest)   = collected weight in viewing directions
c      testd(ntest)   = stdev of weight in viewing directions
c      testn(ntest)   = number of scatterings contributing to testw 
c      testm(ntest,10)= weigt due to different scattering orders 1-9
c                       last element - all the rest

c     other output:
c      npartori = number of particles in position-file
c      npart    = number of particles + replicas of those
c                 crossing the edges of the original region
c      tau      = dynamical (geometric) optical depth
c      ndirect  = number of photons passing directly through

c      outfile  = file to store results of both indirect and direct MC
c                 weights and directions of each escaping photon
c                 +weights in test-direction (indirect MC) 
c      epsfile  = file to store results of both indirect and direct MC
c                 weight of escaping photons tabulated over direction
c      wwwfile  = just weights (indirect MC)

c      npot_mean = mean number of particles the photon may intersect 
c      nbox_mean = mean number of replica boxes studied
c      iscat_mean= mean number of scatterings before escape 

c      ifor     = not used anymore
c-------------------------------------------------------------

c-------------------------------------------------------------
c     Variables transferred via commons to/from main program
c     (just some of them listed here)
c-------------------------------------------------------------

c     mcnew_main.inc:
c     iglobe = whether to use Sun or Saturn as light source (see below)
c     abox,rsaturn = Saturnocentric distance, Saturn radius
c     nlat_sat,nfii_sat,model_sat = divisions of globe, Saturn-model
c     exsat,eysat,ezsat = directions of photons from Saturn

c-------------------------------------------------------------
c     Particle field related variables (in commons in mcnew_sub.inc)

c     xx,yy,zz,ss   position and radii
c     xl,yl,sy      semi-width, tangential shift between adjacent regions
c     npart         number of particles
c     zmin,zmax     minimum and maximum extent of particle field

c-------------------------------------------------------------
c     Photon path related variables (via subroutine calls)

c     px,py,pz     current position
c     ex,ey,ez     current direction
c     ipart_now    particle index of recent scattering                 

c-------------------------------------------------------------
c     LIST OF SUBROUTINES:
c-------------------------------------------------------------

c     READS PARTICLE POSITIONS (called once)
c       subroutine part_read(file)

c     CONSTRUCTS SPATIAL PARTITIONING TABLES (mx,my cells; called once)
c       subroutine tab_init

c     MOVES PHOTON INITIAL POSITION TO ACTUAL CALCULATION REGION (called once) 
c       subroutine retbox(x0,y0)

c     CALCULATES ONE SEGMENT OF PHOTON PATH
c     (to the next scattering or escape of the photon,
c      possibly through several image regions, using path_slice)
c       subroutine path_segment(px,py,pz,ex,ey,ez,ipart,ivis)

c     CALCULATES ONE SLICE OF PHOTON PATH-SEGMENT
c     (in a single region: either the path to next scattering or 
c      to the border, called by path_segment)
c       subroutine path_slice(px,py,pz,ex,ey,ez,ipart,tmin,ivis)

c     CHOOSES PARTICLES WHICH A PHOTON MAY INTERSECT ALONG ITS PATH
c     BASED ON TABLES CREATED IN TAB_INIT
c       subroutine tab_choose(px,py,pz,ex,ey,ez,npot,pottab)

c     CALCULATES SCATTERING: new direction of photon
c       subroutine scatter(px,py,pz,ex,ey,ez,ipart)

c     COLLECTS WEIGHTS ON TEST DIRECTIONS (indirect MC)
c       subroutine indirect(px,py,pz,ipart,iscat)

c     p_ -> CALCULATES PHASE FUNCTION
c       subroutine p_lambert(cose,phasef)
c       subroutine p_hg(cose,phasef,g)
c       subroutine p_hg2(cose,phasef,g1,g2,b)
c       subroutine p_euler(cose,phasef,g)

c     e_ -> CALCULATES PHOTON DIRECTION ANGLES AFTER SCATTERING
c       subroutine e_lambert(sine,cose)
c       subroutine e_hg(sine,cose,g)
c       subroutine e_hg2(sine,cose,g1,g2,b)
c       subroutine e_euler(sine,cose,g)

c     FOR POWER_LAW PHASE FUNCTION, REQUIRES INTERPOLATION TABLES
c       subroutine init_callisto
c       subroutine p_callisto(cose,phasef)
c       subroutine e_callisto(sine,cose)

c     CONSTRUCTION OF PHOTON INITIAL DIRECTIONS FOR SATURN SHINE
c	subroutine saturn_illumination(sunlat,sunfii,nphot_ori,nphot_tod,flux_sat)

c     CPU_TIME COUNTER
c       subroutine tim(cpu)

c     + Additional routines from Numerical Recipes


c------------------------------------------------------------------------
c     LIST OF PHASE-FUNCTIONS INCLUDED (choose by method)
c------------------------------------------------------------------------

c     method

c     = 1  - Lambert surface element's scattering law
c     = 11 - Lambert spherical particle, photon continues from part center
c     = 12 - Lambert spherical particle, continues from scattering point 
c            ppar1,ppar2,ppar3 - not used

c     = 3  - Henyey-Greenstein, photon continues from particle center
c     = 4  - Henyey-Greenstein, photon continues from scattering point 
c            ppar1 = g, anisotropy factor

c     = 5  - Euler, photon continues from particle center
c     = 6  - Euler, photon continues from scattering point 
c            ppar1 = g, anisotropy factor

c     = 13  - 2-lobe HG, photon continues from particle center
c     = 14  - 2-lobe HG, photon continues from scattering point 
c            P = B*P1(g1) + (1-B)*P2(g2)
c            ppar1,ppar2,ppar3 = g1,g2,B

c     =7   - Callisto (n=3.092) photon continues from particle center
c     =8   - Callisto (n=3.092) photon continues from scattering point
c            ppar1,ppar2,ppar3 - not used

c     =9   - Power-law with n, photon continues from particle center
c     =10  - Power-law with n, photon continues from scattering point
c            ppar1=n  (7 and 8 are just special cases of 9 and 10)
c---------------------------------------------------------------------------

c     Some miscellaneous variables:

c     W_m corresponds to sum(W_k *s)
c     testd_apu0(ntest) = W_m for a single photon (in common)
c     testd_apu1(ntest) = collects sum of W_m over photons
c     testd_apu2(ntest) = collects sum of W_m^2 over photons

c     write ss-error to testd
c     testd_apu1_ss(ntest) = collects sum of W_m over photons, ss-only
c     testd_apu2_ss(ntest) = collects sum of W_m^2 over photons, ss-only


c--------------------------------------------------------------

      subroutine mcnew_sub(nphot,posfile,parfile,outfile,
     +           sunfii,sunlat,
     +           ntest0,testfii,testlat,testw0,testd0,testn0,testm0,
     +           idum0,icheck0,mx0,my0,
     +           npart0,npartori0,npot_mean,nbox_mean,iscat_mean,
     +           method0,epsfile,wwwfile,ifor,sredu0,tau0,ndirect,
     +           zredu0,ppar10,ppar20,ppar30)

      include 'mcnew_sub.inc'
      include 'mcnew_main.inc'

      character*60 posfile,parfile,outfile,epsfile,wwwfile
      double precision latsun,fiisun,inidist,lattest,fiitest
      double precision testfii(nnt),testlat(nnt)
      double precision testw0(nnt),testd0(nnt),testn0(nnt)
      double precision testm0(nnt,10)
      double precision iscat_mean,npot_mean,nbox_mean,mu0
      real fiilat(180,90,10)
      double precision testd_corr(nnt),testd_apu1(nnt),testd_apu2(nnt)
c     testd_apu0(nnt) in common

      double precision testd_apu1_ss(nnt),testd_apu2_ss(nnt)

      real ran1

c==============================================================
c       INITIALIZATION
c==============================================================

c==============================================================
c     read auxillary parameter values from parfile
c     parfile='junk' -> use defaults defined here:

      ssca=1.
      maxscat=100
      nbox=10000
      nboxvis=10000

      if(parfile.ne.'junk') then
      open(unit=20,file=parfile,status='old')         
      read(20,*) ssca
      read(20,*) maxscat
      read(20,*) nbox,nboxvis
      close(20)
      endif
      if(nboxvis.lt.0) nboxvis=nbox

c     if nboxvis negative -> force nboxvis equal to nbox
c     nboxvis=0  for testing purposes
c     nboxvis=1  testing: corresponds to mc_atsi3_f.pro
c==============================================================

c     copy variables in call
c     (we want to keep variables both in the call and in the common)
c     (Why? For some reason that has been forgotton!)

      ntest=ntest0
      idum=idum0
      icheck=icheck0
      method=method0
      sredu=sredu0
      zredu=zredu0

      ppar1=ppar10
      ppar2=ppar20
      ppar3=ppar30


c--------------------------------------
c     initialize tables for indirect MC

      do 10 i=1,ntest
         fiitest=ar*testfii(i)
         lattest=ar*testlat(i)
         extest(i)=cos(lattest)*cos(fiitest)
         eytest(i)=cos(lattest)*sin(fiitest)
         eztest(i)=sin(lattest)
         testw(i)=0.d0
         testd(i)=0.d0
         testn(i)=0.d0
         testd_corr(i)=0.d0

         testd_apu1(i)=0.d0
         testd_apu2(i)=0.d0
         testd_apu1_ss(i)=0.d0
         testd_apu2_ss(i)=0.d0

         do 20 j=1,10
            testm(i,j)=0.d0
 20      continue
 10   continue


c------------------------------------------------------------
c     read posdat-file
c     and create duplicate particles for those
c     partially crossing the boundary

      call part_read(posfile)
      tau0=tau

c-----------------------------------------------------
c     initialization for tab-routines
c     divides particles into mx,my cells
c     do not use if mx,my=1,1

      mx=mx0
      my=my0
      if(mx.gt.1.and.my.gt.1) call tab_init

c-----------------------------------------------------
c     miscellaneous counters

      npot_sum=0
      npot_lkm=0
      nbox_sum=0
      nbox_lkm=0
      iscat_sum=0
      npot_mean=0
      nbox_mean=0
      iscat_mean=0
      ndirect=0

c------------------------------------------------------
c     direction of illumination: exsun,eysun,ezsun
c     (ar = degrees_to_radians)

      latsun=ar*sunlat
      fiisun=ar*sunfii
      exsun=cos(latsun)*cos(fiisun)
      eysun=cos(latsun)*sin(fiisun)
      ezsun=sin(latsun)

c------------------------------------------------------
c     open output files:
c     filename='junk' -> do not use

c     outfile:
c     output of weights + all emerging photon directions
      iout=1
      if(outfile.eq.'junk') iout=0

c     epsfile:
c     output of weights + tabulated emerging directions
c     use fixed 2 by 2 degree intervals, 10 orders = fiilat
      ieps=1
      if(epsfile.eq.'junk') ieps=0

      do 5 ifii=1,180
         do 6 ilat=1,90
            do 7 im=1,10
               fiilat(ifii,ilat,im)=0.
 7          continue
 6       continue
 5    continue

c     output of weights to file wwwfile
c     wwwfile now opened by main program


c--------------------------------------------------------
c     initialize power-law phase-functions

      if(method.ge.7.and.method.le.10) then
         power_expo=ppar1
         if(method.eq.7.or.method.eq.8) power_expo=3.092
         call init_callisto
         ppar3=c_n
      endif


         
c==============================================================
c     LOOP OVER PHOTONS
c==============================================================

      call tim(cpu1)
      IF(NPHOT.GE.1) THEN

         nphots=sqrt(1.d0*nphot)
         ixphot=0
         iyphot=0

      do 100 iphot=1,nphot
         
c     w=weight of photon, initialize to 1
c     ipart_now = particle of previous scattering, initialized to undefined

         w=1.d0
         ipart_now=-999

         do 101 i=1,ntest
            testd_apu0(i)=0.
 101     continue


c------------------------------------
c     choose initial photon direction
c        iglobe = 0     -> SUN
c        iglobe = 1,2,3 -> SATURN (ALL,SOUTH,NORTH)
c        in case of SUN: use previously defined exsun,eysun,ezsun
c        in case of SATURN: read from tables exsat,eysat,ezsat
c           constructed in subroutine saturn_illumination, called previously 
c           by the main program

c ---    SATURN?
         if(iglobe.ge.1) then
            exsun=exsat(iphot)
            eysun=eysat(iphot)
            ezsun=ezsat(iphot)
         endif
c------------------------------------
c     choose point where photon path crosses the central plane
c     uniformly distributed in the calculation region
c     with semi-widths xl,yl

         x0=-xl+ran1(idum)*2.d0*xl
         y0=-yl+ran1(idum)*2.d0*yl
         z0=0.d0

c     alternative method, requires nphot=square number
c     divide calculation region into equally spaced cells and shoot
c     one packet/each cell (to use, uncomment cc-lines)

cc         x0=-xl+(ixphot+0.5)*2.d0*xl/nphots
cc         y0=-yl+(iyphot+0.5)*2.d0*yl/nphots

cc         ixphot=ixphot+1
cc         if(ixphot.eq.nphots) then
cc            ixphot=0
cc            iyphot=iyphot+1
cc         endif


c     move initial photon position along its path so that
c     it starts below or above the particle layer

         inidist=10.d0*xl
         if(ezsun.gt.0) inidist=zmax/ezsun*1.01
         if(ezsun.lt.0) inidist=zmin/ezsun*1.01

         px=x0+exsun*inidist
         py=y0+eysun*inidist
         pz=z0+ezsun*inidist
         ex=-exsun
         ey=-eysun
         ez=-ezsun

c     actual starting point in some image-region ->
c     move an image-photon to the actual calculation region

            if(icheck.ge.1) then 
            write(*,*) '-----------------------------'
            write(*,*) 'iphot = ',iphot
            write(*,*),'ex,ey,ez',ex,ey,ez
            write(*,*),'px,py: before retbox',px,py
            endif

         call retbox(px,py)

            if(icheck.ge.1) then 
            write(*,*),'px,py: after  retbox',px,py
            endif

c==============================================================
c     LOOP OVER SCATTERINGS
c==============================================================
         do 200 iscat=1,maxscat  
  
            iscat_sum=iscat_sum+1
            if(icheck.ge.1) then  
               write(*,*) '--> iscat,px,py,pz',iscat,px,py,pz
            endif

c     forward photon to intersection with particle ipart
c     ivis=1  indicates that this is not a visibility test
c     to be performed for indirect MC
         
      ivis=1
      call path_segment(px,py,pz,ex,ey,ez,ipart,ivis)

            if(icheck.ge.1) then 
            write(*,*),'after path_segment: ipart=',ipart
            endif
      
c     ipart < 0 --> escaped (add weight in direct method)
c     scatab: 0 -> direct, 1 -> single , 2,3... -> multiple

c-----------------------
c     photon has escaped?
c     escape without any scattering -> increment ndirect
 
            if(ipart.le.0) then
               if(iscat.eq.1) ndirect=ndirect+1

c     outfile?
               if(iout.eq.1) then
                  lattab(iphot)=asin(ez)
                  fiitab(iphot)=atan2(ey,ex)
                  scatab(iphot)=iscat-1
               endif
c     epsfile?
               if(ieps.eq.1) then
                  if(iscat.gt.1.and.iscat.le.10) then
                     ifii=int((atan2(ey,ex)*ra+180)/2.)+1
                     ilat=int((asin(ez)*ra+90)/2.)+1
                     fiilat(ifii,ilat,iscat)=fiilat(ifii,ilat,iscat)+w
                  endif
               endif
               if(icheck.ge.1) then 
                  write(*,*),'escape: lat,fii,sca ',
     +            lattab(iphot),fiitab(iphot),scatab(iphot)
               endif
               goto 99
            endif

c-----------------------
c     a new scattering particle was found
c     ipart > 0  indicates scattering (add weight in indirect method)

      if(ipart.ge.1) then
      call scatter(px,py,pz,ex,ey,ez,ipart)
      call indirect(px,py,pz,ipart,iscat)
      ipart_now=ipart
      endif

c     collect ss-error
      if(iscat.eq.1) then
         do 103 i=1,ntest
            testd_apu1_ss(i)=testd_apu1_ss(i)+testd_apu0(i)
            testd_apu2_ss(i)=testd_apu2_ss(i)+testd_apu0(i)**2
 103     continue
      endif

c==============================================================
c     end loop over scatterings
 200  continue
c==============================================================
 99   continue

c     collect ms-errors      
         do 102 i=1,ntest
            testd_apu1(i)=testd_apu1(i)+testd_apu0(i)
            testd_apu2(i)=testd_apu2(i)+testd_apu0(i)**2
 102     continue

c==============================================================
c     end loop over photons
 100  continue
c==============================================================

c     nphot gt 0 ?
      ENDIF

c==============================================================
c     calculate tau_n

      NPHOT_APU=1
      TAU_PATH=0.
      IF(NPHOT.GE.1) NPHOT_APU=NPHOT

            call tim(cpu2)
            fdirect=1.d0*ndirect/nphot_APU
            mu0=cos((90.d0-sunlat)*ar)
            IF(NPHOT.GE.1) tau_path=-log(fdirect)
            tau_n=tau_path*mu0

c     speed / 10^3 photons
            cpu_comp=(cpu2-cpu1)/(1.d0*nphot_APU/1000.d0)


c==============================================================
c     OUTPUT RESULTS FOR THIS ILLUMINATION/POS_FILE
c==============================================================

      testw_mean=0.

      if(ntest.ge.1) then

c     mean weight:
         do 500 i=1,ntest
            testw_mean=testw_mean+testw(i)
 500     continue
         testw_mean=testw_mean/ntest


c     error-estimate for total ss+ms-flux - testd_corr:
         do 501 i=1,ntest
            if(testn(i).ge.2) then
               testd_apu1(i)=(testd_apu1(i)/nphot_APU)**2
               testd_apu2(i)=testd_apu2(i)/nphot_APU
               testd_corr(i)=sqrt(nphot*(testd_apu2(i)-testd_apu1(i)))
            endif
 501     continue


c     error-estimate for ss-flux  - testd:
         do 502 i=1,ntest
            if(testn(i).ge.2) then
               testd_apu1_ss(i)=(testd_apu1_ss(i)/nphot_APU)**2
               testd_apu2_ss(i)=testd_apu2_ss(i)/nphot_APU
               apu=nphot*(testd_apu2(i)-testd_apu1(i))
               testd(i)=0.
               if(apu.ge.0) testd(i)=sqrt(apu)
            endif
 502     continue

      endif

      nphot_lkm=nphot
      if(nphot_lkm.eq.0) nphot_lkm=1
      iscat_mean=(1.*iscat_sum/nphot_APU-1.)

      if(npot_lkm.eq.0) npot_lkm=1
      npot_mean=1.*npot_sum/npot_lkm

      if(nbox_lkm.eq.0) nbox_lkm=1
      nbox_mean=1.*nbox_sum/nbox_lkm

      if(icheck.ge.0) then
          write(*,*) 'mean: testw,iscat',testw_mean,iscat_mean
          write(*,*) 'mx,my,npart,npot',mx,my,npart,npot_mean
          write(*,*) 'nbox,nboxvis,nbox_mean',nbox,nboxvis,nbox_mean
       endif
       npart0=npart
       npartori0=npartori

c*******************************************************************
c     to outfile: weights, photon directions
********************************************************************

      if(iout.eq.1) then
      open(unit=20,file=outfile,status='unknown')

C     - method
      versio=1  
C     - no output of wtab
      versio=2
c     - tab3: errors corrected
      versio=3
C     - WRITE(20,*) SIZEMIN,S20,F20
      versio=4
c     - ppar1,ppar2,ppar3
      versio=5

      write(20,*) versio
      write(20,*) posfile
      write(20,*) sunfii,sunlat,nphot
      write(20,*) ssca,maxscat,nbox,nboxvis,idum

      write(20,*) method
      if(versio.ge.5) write(20,*) ppar1,ppar2,ppar3

      write(20,*) ntest
      write(20,*) (testfii(i),i=1,ntest)
      write(20,*) (testlat(i),i=1,ntest)
      write(20,*) (testw(i),i=1,ntest)
      write(20,*) (testn(i),i=1,ntest)
      write(20,*) (testd(i),i=1,ntest)
      write(20,*) (testd_corr(i),i=1,ntest)
      write(20,*) ((testm(i,j),i=1,ntest),j=1,10)

      write(20,*) (fiitab(i),i=1,nphot)
      write(20,*) (lattab(i),i=1,nphot)
      write(20,*) (scatab(i),i=1,nphot)

      write(20,*) iscat_mean,npot_mean,nbox_mean
      write(20,*) npartori,npart,mx,my

c     versio=2
      write(20,*) sredu,zredu,tau
      write(20,*) fdirect,tau_path,tau_n
      write(20,*) cpu2-cpu1,cpu_comp

      WRITE(20,*) SIZEMIN,S20,F20
      close(20)
      endif


c*******************************************************************
c     to epsfile: weights, fiilat-table
********************************************************************

      if(ieps.eq.1) then
      open(unit=20,file=epsfile,status='unknown')

c     versio=1  - no output of wtab
      versio=2
c     tab3: errors corrected
      versio=3
      versio=4
      versio=5

      write(20,*) versio
      write(20,*) posfile
      write(20,*) sunfii,sunlat,nphot
      write(20,*) ssca,maxscat,nbox,nboxvis,idum

c     versio =1 -> method
      write(20,*) method
      if(versio.ge.5) write(20,*) ppar1,ppar2,ppar3

      write(20,*) ntest
      write(20,*) (testfii(i),i=1,ntest)
      write(20,*) (testlat(i),i=1,ntest)
      write(20,*) (testw(i),i=1,ntest)
      write(20,*) (testn(i),i=1,ntest)
      write(20,*) (testd(i),i=1,ntest)
      write(20,*) (testd_corr(i),i=1,ntest)
      write(20,*) ((testm(i,j),i=1,ntest),j=1,10)

      nfii=180
      nlat=90
      nm=10

      write(20,*) nfii,nlat,nm
      write(20,*) fiilat
      write(20,*) iscat_mean,npot_mean,nbox_mean
      write(20,*) npartori,npart,mx,my

c     versio=2
      write(20,*) sredu,zredu,tau
      write(20,*) fdirect,tau_path,tau_n
      write(20,*) cpu2-cpu1,cpu_comp

      WRITE(20,*) SIZEMIN,S20,F20
      close(20)
      endif


c*******************************************************************
c     to wwwfile: weights
********************************************************************

      open(unit=20,file=wwwfile,status='unknown')

c     versio=1  - no output of wtab
      versio=2
c     tab3: errors corrected
      versio=3
      versio=4
      versio=5

      write(20,*) versio
      write(20,*) posfile
      write(20,*) sunfii,sunlat,nphot
      write(20,*) ssca,maxscat,nbox,nboxvis,idum
      write(20,*) method
      if(versio.ge.5) write(20,*) ppar1,ppar2,ppar3

      write(20,*) ntest
      write(20,*) (testfii(i),i=1,ntest)
      write(20,*) (testlat(i),i=1,ntest)
      write(20,*) (testw(i),i=1,ntest)
      write(20,*) (testn(i),i=1,ntest)
      write(20,*) (testd(i),i=1,ntest)
      write(20,*) (testd_corr(i),i=1,ntest)
      write(20,*) ((testm(i,j),i=1,ntest),j=1,10)

      write(20,*) iscat_mean,npot_mean,nbox_mean
      write(20,*) npartori,npart,mx,my

c     versio=2
      write(20,*) sredu,zredu,tau
      write(20,*) fdirect,tau_path,tau_n
      write(20,*) cpu2-cpu1,cpu_comp
      WRITE(20,*) SIZEMIN,S20,F20
      close(20)

c****************************************************************
c     via call: indirect-method     
c****************************************************************

c     note testd_corr is now returned

      do 310 i=1,ntest
         testw0(i)=testw(i)
         testd0(i)=testd_corr(i)
         testn0(i)=testn(i)
         do 320 j=1,10
            testm0(i,j)=testm(i,j)
 320      continue
 310   continue

       return
      end


c**************************************************************

c     SUBROUTINES

c**************************************************************


c*****************************************************
      subroutine part_read(file)
c*****************************************************

c     READS THE POSDAT-FILE
c     ALL POS-VARIABLES TRANSFERRED VIA COMMONS

      include 'mcnew_sub.inc'
      include 'mcnew_main.inc'
      character*60 string1,string2,file

      if(icheck.ge.0) write(*,*) '============================'
      open(unit=10,file=file,status='old')
      read(10,*) string1
      if(icheck.ge.0) write(*,*) 'posdat:',string1
      read(10,*) string2
      read(10,*) orb,xl,yl,sy
      if(icheck.ge.0) write(*,*) 'orb,xl,yl,sy',orb,xl,yl,sy
      read(10,*) npart
      if(icheck.ge.0) write(*,*) 'original npart ',npart

c     xl,yl = radial and tangential semiwidth of calculation region
c     sy = tangential shift between adjacent radial regions
c     make sure |sy| < yl:

      sy_ori=sy
      do 1 i=1,10
         if(sy.gt.yl) sy=sy-2.*yl
         if(sy.lt.-yl) sy=sy+2.*yl
 1    continue
      write(*,*) 'sy: ',sy_ori,sy

      npartori=npart
      lx=2.d0*xl
      ly=2.d0*yl
      zmax=0.d0
      zmin=0.d0

c-------------------------------------------------
c     read particle data
c     optionally: reduce radius, expand z-range

      tau=0.
      do 10 i=1,npart
        read(10,*) ss(i),xx(i),yy(i),zz(i)
        ss(i)=ss(i)*sredu
        zz(i)=zz(i)*zredu
        tau=tau+pii*ss(i)**2
        zzmin=zz(i)-ss(i)
        zzmax=zz(i)+ss(i)
        if(zzmax.gt.zmax) zmax=zzmax
        if(zzmin.lt.zmin) zmin=zzmin
 10   continue
      tau=tau/4./xl/yl
      close(10)

c--------------------------------------------------
c     exclude particles with ss<SIZEMIN?

      if(sizemin.gt.0) then
         tau_mod=0.
         nok=0
         do 15 i=1,npart
            if(ss(i).ge.sizemin) then
               nok=nok+1
               ss(nok)=ss(i)
               xx(nok)=xx(i)
               yy(nok)=yy(i)
               zz(nok)=zz(i)
               tau_mod=tau_mod+pii*ss(nok)**2
            endif
 15      continue
         tau_mod=tau_mod/4./xl/yl

         write(*,*) 'sizemin: TAU, TAU_MOD ',sizemin,tau,tau_mod
         write(*,*) 'sizemin: N, N_MOD ',sizemin,npart,nok
         npart=nok
         tau=tau_mod
      endif

c--------------------------------------------------
c     enhance radius of particles with ss<s20?

      if(s20.gt.0) then
         tau_mod=0.
         do 16 i=1,npart
            if(ss(i).le.s20) then
               ss(i)=ss(i)*f20
            endif
            tau_mod=tau_mod+pii*ss(i)**2
 16      continue
         tau_mod=tau_mod/4./xl/yl
         write(*,*) 's20,f20: TAU, TAU_MOD ',s20,f20,tau,tau_mod
         tau=tau_mod
      endif

c--------------------------------------------------
c     create dublicates of particles whose centers are 
c     close to the borders of the calculation region
c     (these particles extend partially over boundaries, so they may
c     intersect photon paths even via their image particles)
c     These dublicates are treted as normal particles from hereon

        lisa=0
        do 20 i=1,npart
           do 21 in=-1,1
              do 22 im=-1,1
                 if(im.ne.0.or.in.ne.0) then
                    xxt=xx(i)+in*lx
                    yyt=yy(i)+in*sy+im*ly
                    if(abs(xxt).le.xl+ss(i)) then
                       if(abs(yyt).le.yl+ss(i)) then
                          lisa=lisa+1
                          xx(npart+lisa)=xxt
                          yy(npart+lisa)=yyt
                          zz(npart+lisa)=zz(i)
                          ss(npart+lisa)=ss(i)
                       endif
                    endif
                 endif
 22           continue
 21        continue
 20     continue
        npart=npart+lisa
        
        do 30 i=1,npart
           ss2(i)=ss(i)**2
 30     continue
        
        if(icheck.ge.0) write(*,*) 'duplicate edges ',npart
        if(icheck.ge.0) write(*,*) 'zmin,zmax',zmin,zmax
        if(icheck.ge.0) write(*,*) '============================'
        
        return
        end

          
c********************************************************
      subroutine path_segment(px,py,pz,ex,ey,ez,ipart,ivis)
c********************************************************

      include 'mcnew_sub.inc'

c     Calculates the photon path-segment from the present scattering 
c     to the next, which may require following the path through several 
c     image regions (= several "path-slices").
c     The slice of the photon path in a single region is calculated by
c     the subroutine path_slice, which tells whether to advance the path 
c     to the next particle of scattering, or to the boundary of the region.
c     Which is the case is defined by the value of ipart returned by path_slice:
c       ipart=-4,-3,-2,-1    advance to boundary
c       ipart => 1           advance to next particle of scattering = ipart  

c     path-segment calls path_slice repeatedly until either
c      1) next intersection is found,
c      2) photon escapes (z < zmin or z > zmax),
c      3) max number of image searches is reached (ibox=nbox -> escape)

c     arguments of the call of path_segment:
c      px,py,pz = initial photon position -> replaced by the new position  
c      ex,ey,ez = initial photon direction   (not affected in exit)
c      ipart returns the next particle of scattering
c      if ipart <=0 -> photon has escaped

c     Subroutine path_segment is also called by the subroutine indirect,
c     for checking whether the viewing direction is unblocked
c     This in indicated by IVIS=2 (normally IVIS=1)
c      In this case: if ipart => 1, then the viewing direction is blocked

c     nboxtod = number of path-slices (image boxes) checked
c      IVIS=1: use nboxtod=nbox
c      IVIS=2: use nboxtod=nboxvis (different form nbox just for test purposes)  

c     counters (for IVIS=1):
c      nbox_sum = cumulative number of image regions checked 
c      nbox_lkm = number of calls to this program
c      --> used for getting the mean number of image regions checked
 
c---------------------------------------------------------------------------
        ipart=0
        nboxtod=nbox
        if(ivis.eq.2) nboxtod=nboxvis
        if(ivis.eq.1) nbox_lkm=nbox_lkm+1
        if(nboxtod.eq.0) goto 20

c---------------------------------------------------------
c loop over image-boxes

	do 10 ibox=1,nboxtod
	call path_slice(px,py,pz,ex,ey,ez,ipart,tmin,ivis)
        if(ivis.eq.1) nbox_sum=nbox_sum+1

c     ipart_now indicates the previous particle of scattering:
c     on first part_slice (within actual region) use ipart_now
c     defined in the calling routine
c     on next slices (image regions) set to undefined
        ipart_now=-999

c     tmin is the distance along the present path
c     update photon position:

	px=px+ex*tmin
	py=py+ey*tmin
	pz=pz+ez*tmin

        if(icheck.ge.1) then
           if(ivis.eq.1) write(*,*) 'seg: px,py,pz',px,py,pz
        endif

c     Exit from path-segment if an intersection with a particle is found:
	if(ipart.ge.1) goto 20

c     Also exit if the photon is outside the particle layer:
	if(pz.gt.zmax.and.ez.ge.0) goto 20
	if(pz.lt.zmin.and.ez.le.0) goto 20

c     no intersection is found and the photon is still within particle layer:
c     --> move the photon to the next image region
c     This is handled by moving an image-photon to the original region
c      ipart=-4,-3,-2,-1 indicate the boundary which is crossed:
c      each requires are different shift

	if(ipart.eq.-1) then
           px=px-lx
           py=py-sy
	endif
	if(ipart.eq.-2) then
           px=px+lx
           py=py+sy
	endif
	if(ipart.eq.-3) then
           py=py-ly
	endif
	if(ipart.eq.-4) then
           py=py+ly
	endif
	if(py.gt.yl) py=py-ly
	if(py.lt.-yl) py=py+ly

c	continue with new image-box
 10	continue

c	exit from path_segment
 20	continue

	return
	end

c*************************************************************
      subroutine path_slice(px,py,pz,ex,ey,ez,ipart,tmin,ivis)
c*************************************************************

      include 'mcnew_sub.inc'
      dimension tr(4)
      integer pottab(nn)

c     CALCULATES ONE PATH-SLICE = photon path in a single region 
c     (intersection with a particle or with the xy-borders)

c     px,py,pz = photon position   = P
c     ex,ey,ez = photon direction  = E
c     xx,yy,zz = particle position = R

c     ipart    = index of the next particle on the path
c              = or the edge of the region that will be hit next   
c                -1 -> positive x
c                -2 -> negative x
c                -3 -> positive y
c                -4 -> negative y

c----------------------------------------------
c     Search for particles on the current path
c     D = [E * (P-R)]^2 - (P-R)^2 + ss^2  
c     Distance to the particle t=-E*(P-R) - sqrt(D) 
c     intersects if t is real and positive
c     search the particle yielding the smallest positive t

c     particle ipart_now is excluded from the search
c     = particle of previous scattering

      tmin=1.d6
      ipart=0

c-----------------------------------------------------------
c     This same routine is used for checking whether a given
c     direction is blocked by any particle (IVIS=2)
c     In this case exit immediately if any blocking particle
c     is found
c------------------------------------------------------------

c------------------------------------------------------------
c     Search for intersections:

c------------------------------------------------------
c     mx=my=1  -> do not use spatial partition tables ?
c                 all particles are checked
c------------------------------------------------------

      if(mx.le.1.and.my.le.1) then
         do 10 i=1,npart
            if(i.eq.ipart_now) goto 10
            dx=px-xx(i)
            dy=py-yy(i)
            dz=pz-zz(i)
            apu1=-(ex*dx+ey*dy+ez*dz)
            if(apu1.lt.0) goto 10
            
            apu2=dx*dx+dy*dy+dz*dz
            disk=apu1*apu1-apu2+ss2(i)
            if(disk.lt.0.d0) goto 10

            t=apu1-sqrt(disk)
            if(t.le.0) goto 10

            if(ivis.eq.2) then
               ipart=i
               return
            endif

            if(t.lt.tmin) then
               ipart=i
               tmin=t
            endif
 10      continue
      endif

c------------------------------------------------------
c     mx and my > 1  -> use spatial partition tables ?
c     only particles listed in pottab are checked
c------------------------------------------------------

      if(mx.gt.1.and.my.gt.1) then
         call tab_choose(px,py,pz,ex,ey,ez,npot,pottab)
         if(npot.gt.0) then
            do 11 ii=1,npot
               i=pottab(ii)
               if(i.eq.ipart_now) goto 11
               dx=px-xx(i)
               dy=py-yy(i)
               dz=pz-zz(i)
               apu1=-(ex*dx+ey*dy+ez*dz)
               if(apu1.lt.0.d0) goto 11
               
               apu2=dx*dx+dy*dy+dz*dz
               disk=apu1*apu1-apu2+ss2(i)
               if(disk.lt.0.d0) goto 11
               
               t=apu1-sqrt(disk)
               if(t.le.0.d0) goto 11
               
               if(ivis.eq.2) then
                  ipart=i
                  return
               endif

               if(t.lt.tmin) then
                  ipart=i
                  tmin=t
               endif
 11         continue
         endif
      endif


      if(icheck.ge.1) then 
         if(ivis.eq.1) then
            if(ipart.ne.0) write(*,*) 'slice_hit: ipart,tmin',ipart,tmin
         endif
      endif


c     no intersection found ?
c     then search the next border along the photon path

      if(ipart.eq.0) then
         tr(1)=1d6
         tr(2)=1d6
         tr(3)=1d6
         tr(4)=1d6

         if(ex.gt.0) tr(1)=(xl-px)/ex
         if(ex.lt.0) tr(2)=(-xl-px)/ex
         if(ey.gt.0) tr(3)=(yl-py)/ey
         if(ey.lt.0) tr(4)=(-yl-py)/ey

         tmin=2d6
         do 20 i=1,4
            if(tr(i).lt.tmin) then
               ipart=-i
               tmin=tr(i)
            endif
 20      continue
         if(icheck.ge.1) then 
            if(ivis.eq.1) then
               write(*,*) 'slice_edge: ipart,tmin',ipart,tmin
            endif
         endif
      endif

      return
      end


c********************************************************
      subroutine scatter(px,py,pz,ex,ey,ez,ipart)
c********************************************************

      include 'mcnew_sub.inc'
      real ran1

c     Calculate scattering:
c     - update photon direction ex,ey,ez
c     - reduce weight
        w=w*ssca
c     - store current particle index
        ipart_now=ipart

c-------------------------------------------
c     SURFACE-ELEMENT'S LAW ->
c     new direction with respect to unit normal in scattering point
c     scattering point intact

c     SPHERICAL SCATTERING ->
c     new direction with respect to unit normal in incoming direction
c     either keep scattering point intact or
c     move scattering point to particle center,
c     depending on method

c-------------------------------------------------------------------------
c     LAMBERT LAW ?
c     SURFACE-ELEMENT
c-------------------------------------------------------------------------
      if(method.eq.1) then
c------------------------------------
c     unit normal in scattering point
         ax=px-xx(ipart)
         ay=py-yy(ipart)
         az=pz-zz(ipart)
         aa=sqrt(ax**2+ay**2+az**2)
         ax=ax/aa
         ay=ay/aa
         az=az/aa
         sine=sqrt(ran1(idum))
         cose=sqrt(1.d0-sine**2)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     LAMBERT LAW ?
c     SPHERICAL PARTICLES
c-------------------------------------------------------------------------
      if(method.eq.11) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     move scattering point to particle center
         ax=ex
         ay=ey
         az=ez
         px=xx(ipart)
         py=yy(ipart)
         pz=zz(ipart)
         call e_lambert(sine,cose)
         fii=ran1(idum)*2.*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif
c------------------------------------
      if(method.eq.12) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     scattering point intact
         ax=ex
         ay=ey
         az=ez
         call e_lambert(sine,cose)
         fii=ran1(idum)*2.*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     HG LAW ?
c     SPHERICAL PARTICLES
c-------------------------------------------------------------------------
      if(method.eq.3) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     move scattering point to particle center
         ax=ex
         ay=ey
         az=ez
         px=xx(ipart)
         py=yy(ipart)
         pz=zz(ipart)
c     random direction according to hg-law
c     ppar1=g-factor
         call e_hg(sine,cose,ppar1)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif
c------------------------------------
      if(method.eq.4) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     scattering point unaffected
         ax=ex
         ay=ey
         az=ez
         call e_hg(sine,cose,ppar1)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     HG2 LAW ?
c     SPHERICAL PARTICLES
c-------------------------------------------------------------------------
      if(method.eq.13) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     move scattering point to particle center
         ax=ex
         ay=ey
         az=ez
         px=xx(ipart)
         py=yy(ipart)
         pz=zz(ipart)
c     random direction according to hg2-law
c     ppar1=g1-factor
c     ppar2=g2-factor
c     ppar3= b-factor
         call e_hg2(sine,cose,ppar1,ppar2,ppar3)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif
c------------------------------------
      if(method.eq.14) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     scattering point unaffected
         ax=ex
         ay=ey
         az=ez
         call e_hg2(sine,cose,ppar1,ppar2,ppar3)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     EULER LAW ?
c     SPHERICAL PARTICLES
c-------------------------------------------------------------------------
      if(method.eq.5) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     move scattering point to particle center
         ax=ex
         ay=ey
         az=ez
         px=xx(ipart)
         py=yy(ipart)
         pz=zz(ipart)
c     random direction according to euler-law
c     ppar1=x-factor
         call e_euler(sine,cose,ppar1)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif
c------------------------------------
      if(method.eq.6) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     scattering point unaffected
         ax=ex
         ay=ey
         az=ez
         call e_euler(sine,cose,ppar1)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     DONES POWER LAW ?
c     SPHERICAL PARTICLES
c-------------------------------------------------------------------------
      if(method.eq.7.or.method.eq.9) then
c------------------------------------------
c     unit normal in the direction of incoming photon
c     move scattering point to particle center
         ax=ex
         ay=ey
         az=ez
         px=xx(ipart)
         py=yy(ipart)
         pz=zz(ipart)
c     random direction according to
c     DONES et al 1993  power-law
c     n=ppar1
c     iteration tables created in SUBROUTINE power_init
         call e_callisto(sine,cose)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c------------------------------------
      if(method.eq.8.or.method.eq.10) then
c------------------------------------
c     unit normal in the direction of incoming photon
c     scattering point unaffected
         ax=ex
         ay=ey
         az=ez
         call e_callisto(sine,cose)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c-------------------------------------------------------------------------
c     ISOTROPIC (FULL SPHERE)?
      if(method.eq.2) then 
c-------------------------------------------------------------------------
c     unit normal in scattering point
         ax=px-xx(ipart)
         ay=py-yy(ipart)
         az=pz-zz(ipart)
         aa=sqrt(ax**2+ay**2+az**2)
         ax=ax/aa
         ay=ay/aa
         az=az/aa
         cose=ran1(idum)*2.d0-1.d0
         sine=sqrt(1.d0-cose**2)
         fii=ran1(idum)*2.d0*pii
         cosf=cos(fii)
         sinf=sin(fii)
      endif

c---------------------------------------------------------
c     store old direction and
c     calculate new direction

      exold=ex
      eyold=ey
      ezold=ez

      if(cose.lt.-1.d0.or.cose.gt.1.d0) then
         sine=0.d0
         write(*,*) 'error: scatter: cose=',cose
         if(cose.lt.-1.) cose=-.99999999d0
         if(cose.gt. 1.) cose= .99999999d0
      endif

      if(abs(az).ne.1.) then
         apu1=1.d0/(1.d0-az**2)
         apu2=sqrt(1.d0/apu1)
         az2=az*cose+sine*cosf*apu2
         ax2=apu1*(ax*(cose-az*az2)-sine*sinf*ay*apu2)
         ay2=apu1*(ay*(cose-az*az2)+sine*sinf*ax*apu2)
         ex=ax2
         ey=ay2
         ez=az2
      endif
         
      if(abs(az).eq.1.) then
         ez=az/abs(az)*cose
         ex=sine*cosf
         ey=sine*sinf
      endif

      if(icheck.ge.1) then 
         write(*,*) 'scatter: cose,sine',cose,sine   
         write(*,*) 'scatter: OLD',exold,eyold,ezold
         write(*,*) 'scatter: NEW',ex,ey,ez
      endif
 
      return
      end

c********************************************************
      subroutine indirect(px,py,pz,ipart,iscat)
c********************************************************

c     test whether the scattering point is visible from
c     the viewing direction 
c     and if visible add weight to wieving direction

      include 'mcnew_sub.inc'

c-----------------------------------------
c     reference direction for scattering
c-----------------------------------------

c     unit normal at scattering point
      if(method.le.2) then 
         ax=px-xx(ipart)
         ay=py-yy(ipart)
         az=pz-zz(ipart)
         aa=sqrt(ax**2+ay**2+az**2)
         ax=ax/aa
         ay=ay/aa
         az=az/aa
      endif

c     incoming direction
      if(method.ge.3) then
         ax=exold
         ay=eyold
         az=ezold
      endif

c--------------------------------------------------------
c     Check weight in test directions
c--------------------------------------------------------

      do 10 itest=1,ntest
         ipart_now=ipart

         ex2=extest(itest)
         ey2=eytest(itest)
         ez2=eztest(itest)
         cose=ex2*ax+ey2*ay+ez2*az

         if(cose.lt.-1.d0.or.cose.gt.1.d0) then
            write(*,*) 'error: indirect: cose=',cose
            if(cose.lt.-1.d0) cose=-.99999999d0
            if(cose.gt. 1.d0) cose= .99999999d0
         endif


c 1) -----------------------------------------------------
c     scattering point behind the particle of scattering ?
c     only for lambert-surface-element
         if(method.le.1) then
            if(cose.lt.0) goto 10
         endif

c 2) -----------------------------------------------------
c     scattering point blocked by some other particles?
c     ishadow < 1 --> no particles on the way
c     ivis=2   indicates that this is visibility test
c     ex2,ey2,ez2 dummies

         ivis=2
         px2=px
         py2=py
         pz2=pz

         call path_segment(px2,py2,pz2,ex2,ey2,ez2,ishadow,ivis)
         if(ishadow.ge.1) goto 10

c 3) -----------------------------------------------------
c    scattering point is visible from test direction
c     add weight in viewing direction
c     Note: the factor of pi is taken into account in idl-routines
c     -> here phase = real phase-function * pi

c        Lambert-surface-element-phase-function
         if(method.eq.1) phasef=cose

c        Isotropic-phase-function
         if(method.eq.2) phasef=0.25d0

c        Lambert-sphere-phase-function
         if(method.eq.11.or.method.eq.12) then
            call p_lambert(cose,phasef)
         endif

c        HG-sphere-phase-function
c        ppar1=g-factor
         if(method.eq.3.or.method.eq.4) then
            call p_hg(cose,phasef,ppar1)
         endif

c        HG2-sphere-phase-function
c        ppar1=g1-factor
c        ppar2=g2-factor
c        ppar3=B-factor
         if(method.eq.13.or.method.eq.14) then
            call p_hg2(cose,phasef,ppar1,ppar2,ppar3)
         endif

c        EULER-sphere-phase-function
c        ppar1=x-factor
         if(method.eq.5.or.method.eq.6) then
            call p_euler(cose,phasef,ppar1)
         endif

c        POWER-sphere-phase-function
         if(method.ge.7.and.method.le.10) then
            call p_callisto(cose,phasef)
         endif

         if(icheck.ge.1) then 
            write(*,*) 'indirect: cose,phasef',cose,phasef   
         endif

c     weights and errors in indirect method

         testw(itest)=testw(itest)+w*phasef
         testd(itest)=testd(itest)+(w*phasef)**2
         testn(itest)=testn(itest)+1.
         if(iscat.le.9) testm(itest,iscat)=testm(itest,iscat)+w*phasef
         if(iscat.gt.9) testm(itest,10)=testm(itest,10)+w*phasef

c     weight due to a single photon
         testd_apu0(itest)=testd_apu0(itest)+w*phasef

 10   continue
      return
      end

c********************************************************
      subroutine retbox(x0,y0)
c********************************************************

c     moves position px,py into the actual box
c     image points:
c     x -> x + in * LX 
c     y -> y + in * SY + im * LY
c     note the opposite sign as compared to e.g. mosa5_f.pro
c     LX=2*xl
c     LY=2*yl

      include 'mcnew_sub.inc'

      if(x0.ge.0) in=int((x0+xl)/lx)
      if(x0.lt.0) in=-int((-x0+xl)/lx)

      apuy=y0-in*sy
      if(apuy.ge.0) im=int((apuy+yl)/ly)
      if(apuy.lt.0) im=-int((-apuy+yl)/ly)

      x1=x0-in*lx
      y1=y0-in*sy-im*ly
      x0=x1
      y0=y1

      if(abs(x0).gt.xl) write(*,*) 'abs(x0) > xl',x0
      if(abs(y0).gt.yl) write(*,*) 'abs(y0) > yl',y0

      return
      end

c***********************************
      subroutine tab_init
c***********************************

      include 'mcnew_sub.inc'
      dimension ix(nn),iy(nn)

c     construct tables used for limiting the number of
c     particles checked for intersections:
c      divides particles into mx*my subcells
c      cxtab(ic),cytab(ic) = centers of cell ic
c      dtab(ic) = maximum extent of particles in cell ic
c      dtab2=dtab**2
c      indtab contain indices of all particles, 
c      sorted according to subcells:
c      indtab(mintab(ic):mintab(ic)+numtab(ic)) = particles in cell ic

c----------------------------------------------
c     1) search subcell-index for each particle
c     there are particles whose centers are outside the actual region
c     (= dublicate particles in pos_read)
c     -> sub-cells must cover also these particles
c        (this is why for example xmax is not equal to xl)

      xmin=1d6
      ymin=1d6
      xmax=-1d6
      ymax=-1d6
      do 10 i=1,npart
      if(xx(i).gt.xmax) xmax=xx(i)
      if(xx(i).lt.xmin) xmin=xx(i)
      if(yy(i).gt.ymax) ymax=yy(i)
      if(yy(i).lt.ymin) ymin=yy(i)
 10   continue

      xmax=xmax*1.01d0
      xmin=xmin*1.01d0
      ymax=ymax*1.01d0
      ymin=ymin*1.01d0
      dxcell=(xmax-xmin)/mx
      dycell=(ymax-ymin)/my

      do 20 i=1,npart
         ix(i)=int((xx(i)-xmin)/dxcell)+1
         iy(i)=int((yy(i)-ymin)/dycell)+1
 20   continue

c----------------------------------------------
c     2) place indices etc.

      lask=0
      indmui=0
      do 30 i=1,mx
         do 31 j=1,my
            lask=lask+1
            cxtab(lask)=xmin+(i-0.5)*dxcell
            cytab(lask)=ymin+(j-0.5)*dycell
            nummui=0
            ssmui=0.
            do 32 k=1,npart
               if(ix(k).eq.i.and.iy(k).eq.j) then
                  nummui=nummui+1
                  indmui=indmui+1
                  indtab(indmui)=k
                  if(ss(k).gt.ssmui) ssmui=ss(k)
               endif
 32         continue
            numtab(lask)=nummui
            mintab(lask)=indmui-nummui+1
            dtab2(lask)=(dxcell/2.+ssmui)**2+(dycell/2.+ssmui)**2
            dtab(lask)=sqrt(dtab2(lask))
 31      continue
 30   continue

      return
      end

c*********************************************************
      subroutine tab_choose(px,py,pz,ex,ey,ez,npot,pottab)
c**********************************************************

c     choose the subcells where the current photon path may hit a particle
c     go through subcells:
c     - reject the subcell if its center is too far from the photon path
c     - add particles of the subcell to pottab if a potential hit possible
c     in case of small |ez| additional tests

c     returns:
c     pottab= list of particles to be checked
c     npot  = number of particles in pottab 

      include 'mcnew_sub.inc'
      integer pottab(nn)
      integer pottabz(nn)

      npot=0

c---------------------------------------
c     additional test for small |ez|
c      Construct the upper and lower vertical limits 
c      for the current photon path: ztmin, ztmax
c      If these are outside the particle layer defined
c      by zmin,zmax -> npot = 0

c     |ez| limit taken as 0.02 (not necessarily optimal value!)
c     speeds calculations by a factor of about 2 for 1995 RPX.

      if(abs(ez).lt.0.02) then
         dzmax=(xl+yl)*2.*ez
         if(ez.gt.0) then
            ztmax=pz+dzmax
            ztmin=pz
         endif
         if(ez.le.0) then
            ztmin=pz+dzmax
            ztmax=pz
         endif
         if(ztmin.gt.zmax.or.ztmax.lt.zmin) return
      endif

c---------------------------------------
c     calculate projection of path to xy-plane

      ee=sqrt(ex**2+ey**2)
      ex0=ex/ee
      ey0=ey/ee

      do 10 i=1,mx*my
         if(numtab(i).le.0) goto 10
         dx=cxtab(i)-px
         dy=cytab(i)-py
         epd=ex0*dx+ey0*dy
         d2=dx**2+dy**2
         tmin=epd

c     potential hit possible ?
c     check distance in xy-plane
         if(tmin.gt.0) smin2=d2-epd**2   
         if(tmin.le.0) smin2=d2
         if(smin2.gt.dtab2(i)) goto 10

c     check vertical distance 
         if(tmin.gt.dtab(i)) then
            zpot=pz+(tmin-dtab(i))*ez
            if(zpot.gt.zmax.or.zpot.lt.zmin) goto 10
         endif

c     add particles of this subcell to pottab
         do 20 j=1,numtab(i)
            npot=npot+1
            k=mintab(i)+(j-1)
            pottab(npot)=indtab(k)
 20      continue
 10      continue


c---------------------------------------
c     additional test for small |ez|
c     if ztmax-ztmin << zmax-zmin it is advantageous to discard
c     from pottab all particles which fall outside [ztmin,ztmax],
c     even if extra tests are required here

      if(abs(ez).lt.0.02) then
         npotz=npot
         do 30 i=1,npotz
            pottabz(i)=pottab(i)
 30      continue
         npot=0
         do 40 i=1,npotz
            ii=pottabz(i)
            if(zz(ii)-ss(ii).gt.ztmax) goto 36
            if(zz(ii)+ss(ii).lt.ztmin) goto 36
            npot=npot+1
            pottab(npot)=ii
 36         continue
 40      continue
      endif
c---------------------------------------

         npot_sum=npot_sum+npot
         npot_lkm=npot_lkm+1
         return
         end

c****************************************************
      subroutine p_lambert(cose,phasef)
c****************************************************
c     Lambert-sphere: phase-function (multiplied by pi)
      implicit double precision (a-h,o-z)
      pii=3.141592654d0
      vak=2.d0/pii/3.d0

      eps=acos(cose)
      phasef=vak*(sin(eps)-eps*cose)
      return
      end

c****************************************************
      subroutine e_lambert(sine,cose)
c****************************************************
c     Lambert-sphere: photon-direction
c     inte=1./6./!pi*(-3.*sin(2.*alfa)+4.*alfa+2.*alfa*cos(2.*alfa))
      implicit double precision (a-h,o-z)
      real ran1
      pii=3.141592654d0

      apu=ran1(idum)
      eps1=0.
      eps2=pii
      res1=0.
      res2=1.d0
      eps=0.5*(eps1+eps2)

c     bisektion
      vak=1.d0/6.d0/pii
      do 10 ite=1,20
         res=vak*(-3.d0*sin(2.d0*eps)+4.d0*eps+2.d0*eps*cos(2.d0*eps))
         if(res.gt.apu) then
            eps2=eps
            res2=res
         endif
         if(res.le.apu) then
            eps1=eps
            res1=res
         endif
         eps=0.5d0*(eps1+eps2)
 10   continue

      if(eps.lt.0.or.eps.gt.pii) 
     +     write(*,*) 'error: e_lambert: eps=',eps
      sine=sin(eps)
      cose=cos(eps)
      return
      end

c****************************************************
      subroutine p_hg(cose,phasef,g)
c****************************************************
c     Henyey-Greenstein: phase-function (multiplied by pi)
      implicit double precision (a-h,o-z)

      phasef=0.25d0*(1.d0-g**2)/(1.d0-2.d0*g*cose+g*g)**(1.5d0)
      return
      end

c****************************************************
      subroutine e_hg(sine,cose,g)
c****************************************************      
c     Henyey-Greenstein: photon directions
      implicit double precision (a-h,o-z)
      real ran1

      apu=1.d0-2.d0*ran1(idum)
      apu2=apu*apu
      cose=((apu2-1.d0)*g**3+2.d0*apu*g**2+
     +     (apu2+3.d0)*g+2.*apu)/2.d0/(1.d0+apu*g)**2
      sine=sqrt(1.d0-cose**2)
      return
      end

c****************************************************
      subroutine p_hg2(cose,phasef,g1,g2,B)
c****************************************************
c     2-lobe Henyey-Greenstein: phase-function (multiplied by pi)
      implicit double precision (a-h,o-z)

      phasef1=0.25d0*(1.d0-g1**2)/(1.d0-2.d0*g1*cose+g1**2)**(1.5d0)
      phasef2=0.25d0*(1.d0-g2**2)/(1.d0-2.d0*g2*cose+g2**2)**(1.5d0)
      phasef=b*phasef1+(1.-b)*phasef2
      return
      end

c****************************************************
      subroutine e_hg2(sine,cose,g1,g2,b)
c****************************************************
c     2-lobe Henyey-Greenstein:  photon-direction
c     eps=mu
      implicit double precision (a-h,o-z)
      real ran1
      pii=3.141592654

      apu=ran1(idum)
      eps1=1.d0
      eps2=-1.d0
      res1=0.d0
      res2=1.d0
      eps=0.5d0*(eps1+eps2)

c     bisektion
      b1=b*(1.-g1**2)
      b2=(1.d0-b)*(1.-g2**2)

      do 10 ite=1,20
         apu1=b1/4.d0/g1*((1.+g1**2-2.*g1*eps)**(-0.5)-1./(1.-g1))
         apu2=b2/4.d0/g2*((1.+g2**2-2.*g2*eps)**(-0.5)-1./(1.-g2))
         res=-2.d0*(apu1+apu2)
         if(res.gt.apu) then
            eps2=eps
            res2=res
         endif
         if(res.le.apu) then
            eps1=eps
            res1=res
         endif
         eps=0.5*(eps1+eps2)
 10   continue

      if(eps.lt.-1.or.eps.gt.1) 
     +     write(*,*) 'error: e_hg2: eps=',eps
      cose=eps
      sine=sqrt(1.-eps**2)
      return
      end

c****************************************************
      subroutine p_euler(cose,phasef,g)
c****************************************************
c     EULER: phase-function (multiplied by pi)
      implicit double precision (a-h,o-z)
      
      phasef=0.25d0*(1.d0+g*cose)
      return
      end
      
c****************************************************
      subroutine e_euler(sine,cose,g)
c****************************************************
c     EULER:  photon directions
      implicit double precision (a-h,o-z)
      real ran1

      apu=2.d0-4.d0*ran1(idum)
      cose=(g+apu)/(1.d0+sqrt(1.d0+g*(g+apu)))
      sine=sqrt(1.d0-cose**2)
      return
      end

c****************************************************
      subroutine init_callisto
c****************************************************      
c     POWER-LAW
c     Create iteration tables
c     transferred via common      
c     read ppar1 from here:
      include 'mcnew_sub.inc'
      double precision power_apu(1000)
      save

      dalpha=pii/400.d0
      do 10 i=1,401
         power_alpha(i)=pii*(1.d0*i-1.d0)/400.d0
 10   continue
      do 20 i=1,400
         alpha=pii*(1.d0*i-0.5d0)/400.d0
         power_apu(i)=sin(alpha)*alpha**power_expo
 20   continue
         
      summa=0.5d0*power_apu(1)*dalpha
      do 30 i=2,400
         dsumma=0.5d0*(power_apu(i-1)+power_apu(i))*dalpha
         power_inte(i)=summa+0.5d0*dsumma
         summa=summa+dsumma
 30   continue
         power_inte(1)=0.d0
         power_inte(401)=summa+0.5d0*power_apu(400)*dalpha

c     normalize
         c_n=2.d0/power_inte(401)
         do 40 i=1,401
            power_inte(i)=c_n*power_inte(i)
 40      continue
   

        call spline(power_inte,power_alpha,401,1.d30,1.d30,power_dinte)

c     test
         if(icheck.ge.1) then

            write(*,*) 'INIT_CALLISTO: method',method
            write(*,*) 'INIT_CALLISTO: n, c_n',power_expo,c_n
            write(*,*) 'INIT_CALLISTO: alpha ',
     +           power_alpha(1),power_alpha(401)
            write(*,*) 'INIT_CALLISTO: inte  ',
     +           power_inte(1),power_inte(401)         
            open(unit=60,file='SPLINE_TEST',status='unknown')
            write(60,*) power_expo,c_n
            do 60 i=1,401
              write(60,*) i,power_inte(i),power_alpha(i),power_dinte(i)
 60         continue
            close(60)
         endif

         return
         end

c****************************************************
      subroutine p_callisto(cose,phasef)
c****************************************************
      include 'mcnew_sub.inc'
c     POWER_LAW: phase-function (multiplied by pi)

         alpha=acos(cose)
         phasef=0.25d0*c_n*alpha**power_expo

         return
         end

c****************************************************
      subroutine e_callisto(sine,cose)
c****************************************************
      include 'mcnew_sub.inc'
      real ran1
c     POWER-LAW:  photon directions

c        returns y=alpha vs x=2.*r
         apu=2.d0*ran1(idum)
         call splint(power_inte,power_alpha,power_dinte,401,apu,alpha)

         if(alpha.lt.0.or.alpha.gt.pii) 
     +   write(*,*) 'error: e_callisto: eps=',alpha
         cose=cos(alpha)
         sine=sqrt(1.d0-cose**2)

         return
         end

c***************************************************
	SUBROUTINE SATURN_ILLUMINATION
     +  (sunlat,sunfii,nphot_ori,nphot_tod,flux_sat)
c***************************************************

c     given sunlat,sunfii,nphot,r/a returns tables
c     exsat,eysat,ezsat = direction of incoming photons
c     --> returned via common

c     nphot_ori = goal number of photons
c     nphot_tod = actual number
c     flux_tab(fii,lat) = flux from each surface element

	implicit double precision (a-h,o-z)
	double precision flux_tab(400,200)
	double precision latsun,lat,mu,mu0
	include 'mcnew_main.inc'
	
	dimension xxtab(6)
	dimension atab1(6),atab2(6),atab3(6),atab4(6)
	dimension btab1(6),btab2(6),btab3(6),btab4(6)
	dimension datab1(6),datab2(6),datab3(6),datab4(6)
	dimension dbtab1(6),dbtab2(6),dbtab3(6),dbtab4(6)

	data xxtab/0.,30.,60.,90.,120.,150./

c	BELTS, BLUE
	data atab1/0.6736, 0.6197, 0.5726, 0.5948, 0.7869, 2.1323/
	data btab1/1.1339, 1.1254, 1.1552, 1.1970, 1.2385, 1.4539/
c	ZONES, BLUE
	data atab2/0.5950, 0.5659, 0.5406, 0.5290, 0.6015, 1.5886/
	data btab2/1.0948, 1.0938, 1.1372, 1.1606, 1.1636, 1.3727/
c	BELTS, RED
	data atab3/1.6854, 1.5946, 1.4529, 1.3386, 1.3653, 2.2539/
	data btab3/1.4812, 1.4779, 1.4586, 1.4224, 1.3637, 1.3385/
c	ZONES, RED
	data atab4/1.6911, 1.5995, 1.4567, 1.3436, 1.3680, 2.2021/
	data btab4/1.4825, 1.4794, 1.4592, 1.4221, 1.3618, 1.3341/

	ar=0.017453293d0
	ra=57.29577951d0
	pi=3.141592654d0

	call tim(cpu1)
	itest=0
	if(itest.eq.1) open(unit=50,file='DEBUG',status='unknown')

c***************************************************
c       MODEL_SAT= 1 Lambert-planet 
c       MODEL_SAT=>2  -> need tabulation for 
c       Barkstrom-law (from Dones et al. 1993)

	if(model_sat.gt.1) then
	call spline(xxtab,atab1,6,1d30,1d30,datab1)
	call spline(xxtab,atab2,6,1d30,1d30,datab2)
	call spline(xxtab,atab3,6,1d30,1d30,datab3)
	call spline(xxtab,atab4,6,1d30,1d30,datab4)
	call spline(xxtab,btab1,6,1d30,1d30,dbtab1)
	call spline(xxtab,btab2,6,1d30,1d30,dbtab2)
	call spline(xxtab,btab3,6,1d30,1d30,dbtab3)
	call spline(xxtab,btab4,6,1d30,1d30,dbtab4)
	endif

        if(itest.eq.1) then
           do 5 i=1,181
              alpha=i-1.
              call splint(xxtab,atab1,datab1,6,alpha,aa1)
              call splint(xxtab,atab2,datab2,6,alpha,aa2)
              call splint(xxtab,atab3,datab3,6,alpha,aa3)
              call splint(xxtab,atab4,datab4,6,alpha,aa4)
              write(50,*) alpha,aa1,aa2,aa3,aa4
 5         continue
		 
           write(50,*) ''
           write(50,*) 'alpha,bb1,bb2,bb3,bb4'

           do 6 i=1,181
              alpha=i-1.
              call splint(xxtab,btab1,dbtab1,6,alpha,aa1)
              call splint(xxtab,btab2,dbtab2,6,alpha,aa2)
              call splint(xxtab,btab3,dbtab3,6,alpha,aa3)
              call splint(xxtab,btab4,dbtab4,6,alpha,aa4)
              write(50,*) alpha,aa1,aa2,aa3,aa4
 6         continue
        endif


c***************************************************
c INTEGRATION OF FLUX over fii,lat grid on the planet
c (ignore flattening)
c***************************************************

	nlat=nlat_sat
	nfii=nfii_sat
	dlat=pi/nlat
	dfii=2.d0*pi/nfii

	if(itest.eq.1) then
	   WRITE(50,*) 'NLAT,NFII',NLAT,NFII
	   WRITE(50,*) 'DLAT,DFII',DLAT,DFII
	   WRITE(50,*) 'ABOX,RSATURN,PI',ABOX,RSATURN,PI
	endif

c       direction of sun

	latsun=ar*sunlat
	fiisun=ar*sunfii
	sx=cos(latsun)*cos(fiisun)
	sy=cos(latsun)*sin(fiisun)
	sz=sin(latsun)

	      if(itest.eq.1) then
		 WRITE(50,*) 'SUNLAT,SUNFII',SUNLAT,SUNFII
		 WRITE(50,*) 'SX,SY,SZ',SX,SY,SZ
	      endif

c      integration over Saturn's globe
c       iglobe=1 -> all
c       iglobe=2 -> SOUTH ONLY
c       iglobe=3 -> NORTH ONLY

	FLUX_SAT=0.
	do 10 ilat=1,nlat
	   lat=-pi/2.d0+(ilat-.5d0)*dlat
	   if(iglobe.eq.2.and.lat.ge.0) goto 10
	   if(iglobe.eq.3.and.lat.le.0) goto 10
	   cosl=cos(lat)
	   sinl=sin(lat)

	   do 20 ifii=1,nfii
	      flux_tab(ifii,ilat)=0.
	      fii=(ifii-0.5)*dfii
	      if(itest.eq.1) then
		 WRITE(50,*) 'ILAT,IFII',ILAT,IFII,LAT,FII
	      endif

c       point on the planet
	      px=cosl*cos(fii)
	      py=cosl*sin(fii)
	      pz=sinl

c      direction to SUN
	      mu0=(px*sx+py*sy+pz*sz)
	      if(itest.eq.1) then
		 WRITE(50,*) 'PX,PY,PZ,MU0',PX,PY,PZ,MU0
	      endif
	      if(mu0.le.0) goto 20

c       direction to RING
	      dx=abox-rsaturn*px
	      dy=-rsaturn*py
	      dz=-rsaturn*pz
	      dd=sqrt(dx*dx+dy*dy+dz*dz)
	      mu =(px*dx+py*dy+pz*dz)/dd
	      if(itest.eq.1) then
		 WRITE(50,*) 'DX,DY,DZ,MU',DX,DY,DZ,MU
	      endif
	      if(mu.le.0) goto 20

c       surface element - multiply here with mu
c       and dz/dd = cos of incidence on ring
	      da=cosl*dlat*dfii*rsaturn**2
	      dfac=da*mu/dd/dd*abs(dz)/dd
	      if(itest.eq.1) then
		 WRITE(50,*) 'DA,DFAC',DA,DFAC
	      endif

c----------------------------
c       MODEL FOR SATURN
c             model_sat=1 -> Lambert
c             model_sat=2 -> RED      0.5*zones+0.5*belts
c             model_sat=3 -> BLUE     0.5*zones+0.5*belts
c             model_sat=4 -> RED+BLUE 0.5*zones+0.5*belts
c             model_sat=5 -> ZONES    0.75*blue+0.25*red
c             model_sat=6 -> BELTS    0.75*blue+0.25*red
c             model_sat=11-> BELTS    RED
c             model_sat=12-> BELTS    BLUE
c             model_sat=13-> ZONES    RED
c             model_sat=14-> ZONES    BLUE
c----------------------------
c       LAMBERT-SATURN (ALB=1)
	      if(model_sat.eq.1) then 
		 flux_tab(ifii,ilat)=mu0*dfac
	      endif

              cosa=(sx*dx+sy*dy+sz*dz)/dd
              alpha=acos(cosa)*ra
c----------------------------
c       RED
	      if(model_sat.eq.2) then 
		 call splint(xxtab,atab3,datab3,6,alpha,aa)
		 call splint(xxtab,btab3,dbtab3,6,alpha,bb)
		 apu1 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 call splint(xxtab,atab4,datab4,6,alpha,aa)
		 call splint(xxtab,btab4,dbtab4,6,alpha,bb)
		 apu2 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 flux_tab(ifii,ilat)=(apu1+apu2)/2.*dfac
	      endif
c----------------------------
c       BLUE
	      if(model_sat.eq.3) then 
		 call splint(xxtab,atab1,datab1,6,alpha,aa)
		 call splint(xxtab,btab1,dbtab1,6,alpha,bb)
		 apu1 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 call splint(xxtab,atab2,datab2,6,alpha,aa)
		 call splint(xxtab,btab2,dbtab2,6,alpha,bb)
		 apu2 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 flux_tab(ifii,ilat)=(apu1+apu2)/2.*dfac
	      endif
c----------------------------
c       RED+BLUE
	      if(model_sat.eq.4) then 
		 call splint(xxtab,atab3,datab3,6,alpha,aa)
		 call splint(xxtab,btab3,dbtab3,6,alpha,bb)
		 apu3 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 call splint(xxtab,atab4,datab4,6,alpha,aa)
		 call splint(xxtab,btab4,dbtab4,6,alpha,bb)
		 apu4 = AA/mu *(mu*mu0/(mu+mu0))**BB
		 call splint(xxtab,atab1,datab1,6,alpha,aa)
		 call splint(xxtab,btab1,dbtab1,6,alpha,bb)
		 apu1 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 call splint(xxtab,atab2,datab2,6,alpha,aa)
		 call splint(xxtab,btab2,dbtab2,6,alpha,bb)
		 apu2 = AA/mu *(mu*mu0/(mu+mu0))**BB  
		 flux_tab(ifii,ilat)=(apu1+apu2+apu3+apu4)/4.*dfac
	      endif
c----------------------------
c       0.75*BLUE+0.25*RED: zones
	      if(model_sat.eq.5) then 
c       red-zones
		 call splint(xxtab,atab4,datab4,6,alpha,aa)
		 call splint(xxtab,btab4,dbtab4,6,alpha,bb)
		 apu4 = AA/mu *(mu*mu0/(mu+mu0))**BB
c       blue-zones
		 call splint(xxtab,atab2,datab2,6,alpha,aa)
		 call splint(xxtab,btab2,dbtab2,6,alpha,bb)
		 apu2 = AA/mu *(mu*mu0/(mu+mu0))**BB 
 		 flux_tab(ifii,ilat)=(0.25*apu4+0.75*apu2)*dfac
	      endif
c----------------------------
c       0.75*BLUE+0.25*RED: belts
	      if(model_sat.eq.6) then 
c       red-belts
		 call splint(xxtab,atab3,datab3,6,alpha,aa)
		 call splint(xxtab,btab3,dbtab3,6,alpha,bb)
		 apu4 = AA/mu *(mu*mu0/(mu+mu0))**BB
c       blue-belts
		 call splint(xxtab,atab1,datab1,6,alpha,aa)
		 call splint(xxtab,btab1,dbtab1,6,alpha,bb)
		 apu2 = AA/mu *(mu*mu0/(mu+mu0))**BB 
		 flux_tab(ifii,ilat)=(0.25*apu4+0.75*apu2)*dfac
	      endif
c----------------------------
c       11: belts-red
	      if(model_sat.eq.11) then 
		 cosa=(sx*dx+sy*dy+sz*dz)/dd
		 alpha=acos(cosa)*ra
		 call splint(xxtab,atab3,datab3,6,alpha,aa)
		 call splint(xxtab,btab3,dbtab3,6,alpha,bb)
		 apu = AA/mu *(mu*mu0/(mu+mu0))**BB
		 flux_tab(ifii,ilat)=apu*dfac
	      endif

c----------------------------
c       12: belts-blue
	      if(model_sat.eq.12) then 
		 cosa=(sx*dx+sy*dy+sz*dz)/dd
		 alpha=acos(cosa)*ra
		 call splint(xxtab,atab1,datab1,6,alpha,aa)
		 call splint(xxtab,btab1,dbtab1,6,alpha,bb)
		 apu = AA/mu *(mu*mu0/(mu+mu0))**BB
		 flux_tab(ifii,ilat)=apu*dfac
	      endif
c----------------------------
c       13: zones-red
	      if(model_sat.eq.13) then 
		 cosa=(sx*dx+sy*dy+sz*dz)/dd
		 alpha=acos(cosa)*ra
		 call splint(xxtab,atab4,datab4,6,alpha,aa)
		 call splint(xxtab,btab4,dbtab4,6,alpha,bb)
		 apu = AA/mu *(mu*mu0/(mu+mu0))**BB
		 flux_tab(ifii,ilat)=apu*dfac
	      endif
c----------------------------
c       14: zones-blue
	      if(model_sat.eq.14) then 
		 cosa=(sx*dx+sy*dy+sz*dz)/dd
		 alpha=acos(cosa)*ra
		 call splint(xxtab,atab2,datab2,6,alpha,aa)
		 call splint(xxtab,btab2,dbtab2,6,alpha,bb)
		 apu = AA/mu *(mu*mu0/(mu+mu0))**BB
		 flux_tab(ifii,ilat)=apu*dfac
	      endif
c----------------------------

	      FLUX_SAT=FLUX_SAT+FLUX_TAB(ifii,ilat)

	      if(itest.eq.1) then
		 WRITE(50,*) 'FLUX_TAB',FLUX_TAB(IFII,ILAT)
		 WRITE(50,*) 'FLUX_SAT',FLUX_SAT
	      endif

 20	   continue
 10	continue


c------------------------------------
c       construct photon direction
c------------------------------------

	lask=0
	if(flux_sat.gt.0) then
	do 30 ilat=1,nlat
	   lat=-pi/2.d0+(ilat-.5d0)*dlat
	   cosl=cos(lat)
	   sinl=sin(lat)

	   do 40 ifii=1,nfii
	      fii=(ifii-0.5)*dfii
	      if(flux_tab(ifii,ilat).gt.0) then 
		 napu=nphot_ori*flux_tab(ifii,ilat)/flux_sat

		 if(napu.ge.1) then 
c       direction to RING
		    px=cosl*cos(fii)
		    py=cosl*sin(fii)
		    pz=sinl
		    dx=abox-rsaturn*px
		    dy=-rsaturn*py
		    dz=-rsaturn*pz
		    dd=sqrt(dx*dx+dy*dy+dz*dz)

		    do 50 ii=1,napu
		       lask=lask+1
		       exsat(lask)=-dx/dd
		       eysat(lask)=-dy/dd
		       ezsat(lask)=-dz/DD
c		       WRITE(50,*) 
c     + 'LASK,EX,EY,EZ',LASK,EXSAT(LASK),EYSAT(LASK),EZSAT(LASK)
 50		    continue
		 endif
	      endif
 40	   continue
 30	continue
	endif

	nphot_tod=lask

	call tim(cpu2)
	if(itest.eq.1) then
	   write(50,*) 'SATURN: fii,flux,nphot,cpu',
     +     sunfii,flux_sat,nphot_tod,cpu2-cpu1
	   close(50)
	endif
	write(*,*) 'SATURN: fii,flux,nphot,cpu',
     +  sunfii,flux_sat,nphot_tod,cpu2-cpu1
        
        return
	end


C***********************************
      SUBROUTINE TIM(cpu2)
C***********************************
c     return elapsed cpu-time in seconds

      double precision cpu2
        real tarray(2)
        cpu=etime(tarray)
        cpu2=cpu
      RETURN
      end


c***********************************************************

c       Auxillary routines from NUMERICAL RECIPES


c**************************************************************
        FUNCTION RAN1(IDUM)
c**************************************************************
c       creates uniformly distributed pseudorandom numbers
c       rnd(1)

      save
      DIMENSION R(97)
      PARAMETER (M1=259200,IA1=7141,IC1=54773,RM1=3.8580247E-6)
      PARAMETER (M2=134456,IA2=8121,IC2=28411,RM2=7.4373773E-6)
      PARAMETER (M3=243000,IA3=4561,IC3=51349)
      DATA IFF /0/
      IF (IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IX1=MOD(IC1-IDUM,M1)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX2=MOD(IX1,M2)
        IX1=MOD(IA1*IX1+IC1,M1)
        IX3=MOD(IX1,M3)
        DO 11 J=1,97
          IX1=MOD(IA1*IX1+IC1,M1)
          IX2=MOD(IA2*IX2+IC2,M2)
          R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
11      CONTINUE
        IDUM=1
      ENDIF
      IX1=MOD(IA1*IX1+IC1,M1)
      IX2=MOD(IA2*IX2+IC2,M2)
      IX3=MOD(IA3*IX3+IC3,M3)
      J=1+(97*IX3)/M3
      IF(J.GT.97.OR.J.LT.1) then
      write(*,*) 'ran1:j',j
        j=1
        endif
      RAN1=R(J)
      R(J)=(FLOAT(IX1)+FLOAT(IX2)*RM2)*RM1
      RETURN
      END


c*******************************************************************
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
c*******************************************************************

c     tabulated function y=f(x) with n elements 
c     yp1,ypn = derivatives at the end points
c     y2 = second derivative of f
	
      INTEGER n,NMAX
      double precision yp1,ypn,x(n),y(n),y2(n)
      PARAMETER (NMAX=500)
      INTEGER i,k
      double precision p,qn,sig,un,u(NMAX)
      if (yp1.gt..99e30) then
        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 11 i=2,n-1
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
        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
11    continue
      if (ypn.gt..99e30) then
        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 12 k=n-1,1,-1
        y2(k)=y2(k)*y2(k+1)+u(k)
12    continue
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.

c*******************************************************************
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
c*******************************************************************

c     tabulated function ya=f(xa) with n elements 
c     y2a = second derivative of f
c     y = interpolated value of f(x)

      INTEGER n
      double precision x,y,xa(n),y2a(n),ya(n)
      INTEGER k,khi,klo
      double precision a,b,h
      klo=1
      khi=n
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'
      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.
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software 5.







