cobLateral.f

Go to the documentation of this file.

c     ******************************************************************
c     *                                                                *
c     * mfoblt: observe lateral distribution of a.s generated by mfobas*
c     *                                                                *
c     ******************************************************************
c
c    /usage/       call mfoblt(zs, zobas, lvl)
c
c  -- input --
c         zs:  slant depth of electron position (g/cm**2)
c      zobas:  slatn depth of observation position //
c        lvl:  level index of depth
c  /$mfoo/ age, eno, xa, ya, w1a, w2a
c  /$mfoi/ rsyma2,
c  /$mfci/ hrasi, hraso, szarcr, rcrosm, w3inp
c  /$mfco/ inas1, inas2, rmol2
c  /$mfch/ blat
c
c  -- process --
c    1)  a.s. center is obtained assuming the incident electron keeps
c        current directon upto zobas
c    2)  establish a.s is symmetric around the center of arrays or not.
c    3)  lateral distribution is taken as density map
c    4)  at distances > rcrosm, lateral distribution is taken by arrays
c        distributed like cross.  arrays are equi-distanced in log10
c        scale.
c
c  -- output --
c     /$mfoo/ xas, yas, rmoll
c     /$mfco/ sizei, sizeo, arhorp, arhorm, arverp, arverm
c
c                                =   =   =   =
c
                       subroutine mfoblt(zs, zobas, lvl)
c
c
       include  'Zmfparm.f'
       include  'Zmfoo.f'
       include  'Zmfoi.f'
       include  'Zmfco.f'
       include  'Zmfci.f'
       include  'Zmfch.f'
c
c
      logical symm
c
c
c
c          get distance between z and zobas
      dst = dstair(w3inp, zs, zobas)
c          a.s. core
      xas=xa + dst*w1a
      yas=ya + dst*w2a
c            molier unit at the level (this may be the one 2 r.l above
c            the level if blat contains '2'.
c            and may be halved one if blat contains 'h')
      rmoll=rmol2(lvl)
c            (distance from center)**2
      ras2=xas**2+yas**2
      if(ras2 .lt. rsyma2) then
c           a.s. core is very close to center; regard it symmetric
         symm=.true.
      else
         symm=.false.
      endif
      if(sqary) then
c                      detect a.s. as densitymap
         call mfrgas(sizei(1,1,lvl), hrasi, lszi, symm)
      endif
c
c         lateral by cross:
c
      if(crsary) then
           call mflcra( arhorp(1,lvl),   arhorm(1,lvl), arverp(1,lvl),
     *     arverm(1,lvl), ncrsar,  szarcr, rmoll, rcrosm, eno, symm)
      endif
c
      end
c     ******************************************************************
c     *                                                                *
c     * mfrgas: take lateral distribution of a.s with mesh-like arrays *
c     *                                                                *
c     ******************************************************************
c
c    /usage/     call mfrgas(array,hspan,nm, symm)
c
c
c  -- input --
c     hspan: half of a.s array span (assumed to be rectangular)
c        nm: no. of arrays is nm*nm
c            arrays are numbered from 1 to nm in x- as well as y-
c            direction.
c      symm: logical.  if t, a.s is regarded as symmetric around center
c  /$mfoo/ eno, rmoll
c
c  -- process --
c    1) for each array, compute distance from center of a.s.  compute ne
c       in each array and count it up
c    2) in doing so, if a.s. can be regarded as symmetric around array
c       center, counting-up of ne is done by using symmetry condition
c
c  -- output --
c     array:  2 dimensional a.s array of which (i,j) component
c             represents unit array of size hspan*2/nm.  the center of
c             the array span is assumed to be on the 1ry direction.
c             no. of electrons falling in each unit is counted up.
c             dimension is (nm,nm).
c
c
c
c                               =   =   =   =
c
      subroutine mfrgas(array, hspan, nm, symm)
c
c
c       -inc $mfparm
c       -inc $mfoo
       include  'Zmfparm.f'
       include  'Zmfoo.f'
       include  'Zmfch.f'
c
c
c
c
      dimension    array(nm,nm)
      logical symm
c
c
c
c             half size of unit mesh ...
      harrys=hspan/nm
c             unit array size
      arrys=harrys * 2
c             unit array size in m.u.
      arrysm=arrys/rmoll
c
c             init of y coordinate of mesh ...
      ymsh=hspan+harrys
c             set index for array counter
      if(symm) then
          mxindx=nm/2
          if(mod(nm,2) .ne. 0) mxindx=mxindx+1
      else
          mxindx=nm
      endif
c
c                            count ne in each array
       do   j=1,mxindx
         if(symm) then
            i1=j
         else
            i1=1
         endif
c            array y position
         ymsh=ymsh-arrys
c
c             init of x coordinate of mesh
         xmsh=-hspan+(i1-2)*arrys+harrys
c
             do   i=i1,mxindx
c                  array x position
                xmsh=xmsh+arrys
c
c                      get relative ne in array(i,j):       enr
                call mfnems(harrys, arrysm, xmsh, ymsh, enr)
c                      ne in array
                an=eno*enr
c                      count-up
                array(i,j)=array(i,j)+an
                if(symm) then
c                      use symmetry condition
c
c
c
c               symmetry position and indexes
c
c              *                !                   -->i
c                *     a        !         b    .
c                   * (i,j)     !  (m-i+1,j) .
c                h    *         !         .   c
c                (j,i)   *      !      .
c                           *   !   .   (m-j+1,i)
c             __________________!_________________
c                           .   !  *(m-j+1,m-i+1)
c               (j,m-i+1).      !     *
c                 g    .        !        *   d
c                   . (i,m-j+1) !(m-i+1,m-j+1)
c                .     f        !           e  *
c              .                !                 *
c             !
c             j
c            *** logic for symm case has been tested (82.07.16)
c
                     k1=nm-i+1
                     k2=nm-j+1
c                        avoid overlapped counting
                     if(i .ne. k1) then
c                               /b/
                            array(k1,j)=array(k1,j)+an
c                               /d/
                            array(k2,k1)=array(k2,k1) + an
c                               /f/
                            array(i,k2)=array(i,k2) + an
                     endif
                     if(i .ne. j) then
c                            /c/
                         array(k2,i)=array(k2,i) + an
c                            /g/
                         array(j, k1)=array(j, k1) + an
c                            /e/
                         array(k1,k2)=array(k1,k2) + an
c                            /h/; for /h/ overlap with /g/ also
c                                 avoided
                         if(i .ne. k1)
     *                       array(j,i)=array(j,i) + an
                     endif
                endif
             enddo
       enddo
      end
c     ******************************************************************
c     *                                                                *
c     * mfnems: computes electron no. in a given square mesh using     *
c     *         n-k function.  electron no. is normalized to give 1    *
c     *         when integralated over area.                           *
c     *                                                                *
c     ******************************************************************
c
c    /usage/        call mfnems(ha, am, xmsh, ymsh, enr)
c
c
c  -- input --
c        ha:  half of array size (in cm)
c        am:  array size in m.u.
c      xmsh:  x of the center of the array, measured in 1ry system (cm)
c      ymsh:  y         //
c   /$mfoo/   age, xas, yas, rmoll
c
c  -- process --
c    1) using n-k or greisen's formula or formula by m.c
c       for lateral structure function,
c       normalized electron no. is  computed which falls in a given
c       mesh area.
c
c  -- output --
c       enr:  normalized electron no. falling in the given mesh.
c
c
c
c                               =   =   =   =
c
                 subroutine mfnems(ha, am, xmsh, ymsh, enr)
c
c       -inc $mfparm
c       -inc $mfoo
       include  'Zmfparm.f'
       include  'Zmfoo.f'
       include  'Zmfch.f'
c
      data sqrtpi/1.77246/
c
c
c          distance to a.s. core from the current mesh (m.u.)
c
c
      r=sqrt( (xmsh-xas)**2 + (ymsh-yas)**2 ) /rmoll
      if(abs(xmsh-xas) .gt. ha    .or.  abs(ymsh-yas) .gt. ha) then
c             a.s core is outside of mesh
         if(blat .eq. 'mc') then
             call mnkg(age, r, rho)
         else
             call nkg(age, r, rho)
         endif
         enr=rho*am*am
      else
c            a.s core is inside of mesh
         if(blat .eq. 'mc') then
             call mnkgsf(age, r0, a, b, cs)
             rr=(r+am/sqrtpi)/r0
             enr=cs*2/(age+2.-a)*am*am* rr**(age-a)
         else
             rr=(r+am/sqrtpi)
             cs=csnkg(age)
             enr=cs*2/age* am*am * rr**(age-2.)
         endif
      endif
      end
c     ******************************************************************
c     *                                                                *
c     * mflcra: takes lateral distribution of a.s. with cross-like     *
c     *         arrayes.  each array is separated equally on log10     *
c     *         scale  (1 scale is devided by 5).                      *
c     *                                                                *
c     ******************************************************************
c
c     /usage/   call mflcra(ahp, ahm, avp, avm, nm, aa, rmoll,
c              * r1, eno, symm)
c
c  -- input --
c      nm:  no. of arrays in each wing. (see fig.)
c      aa:  unit array size (in cm)
c      r1:  minimum array position from center of
c           1ry direction (in cm).  must be > 0
c   rmoll:  m.u. at the level
c     eno:  total electron no. of the shower
c symm:  if symmetrical shower, t
c
c                             c
c                             * nm
c
c                             *  2
c                             *  1
c
c             b   :   :  :    +    :  :   :  a
c                        !-r1-!    1  2   nm
c                             *
c                             *
c
c                             *
c                             d
c
c  -- process --
c    1) using mfnems, electron no. is computed in each array and counted
c       up.
c
c  -- output --
c       ahp:    to express arrays in a-direction
c       ahm:     //                  b
c       avp:     //                  c
c       avm:     //                  d
c
c
c                               =   =   =   =
c
        subroutine mflcra(ahp, ahm, avp, avm, nm, aa, rmoll,
     *  r1, eno, symm)
c
c
c
      dimension ahp(nm),ahm(nm),avp(nm),avm(nm)
      logical symm
c
c              10**.2   array position step(equi-distance in log10)
      data dz/1.58489/
c
c         half of array size
      harys=aa/2
c         array size in m.u
      arysm=  aa/rmoll
c         first array position (x or y)
      z = r1
c
       do   k=1,nm
c
c                   array in a-direction
           call mfnems(harys, arysm, z, 0., enr)
           ahp(k)=ahp(k) + enr*eno
c                   array in b-direction
           if(.not. symm) then
               call mfnems(harys, arysm, -z, 0.,  enr)
           endif
           ahm(k)=ahm(k) + enr*eno
c                   array in c-direction
           if(.not. symm) then
               call mfnems(harys, arysm, 0., z,  enr)
           endif
           avp(k)=avp(k) + enr*eno
c                   array in d-direction
           if(.not. symm) then
               call mfnems(harys, arysm, 0., -z, enr)
           endif
           avm(k)=avm(k) + enr*eno
c
           z =z * dz
       enddo
      end