text:
      SUBROUTINE M50MDT(TSEC50,KFLAG,ROTMTX)
      IMPLICIT REAL*8(A-H,O-Z)
C
C  THIS ROUTINE COMPUTES THE MATRICES USED TO ROTATE VECTORS BETWEEN
C  MEAN OF 1950.0 AND MEAN OF DATE COORDINATES.
C
C  MEAN OF 1950.0 IS A COORDINATE SYSTEM ASSOCIATED WITH THE
C  BEGINNING OF THE BESSELIAN YEAR 1950. IT HAS JULIAN EPHEMERIS
C  DATE 2433282.423, 0.077 DAYS BEFORE EPHEMERIS DATE 1/1/50, 0.0 HRS.
C
C  VARIABLE DIM TYPE I/O DESCRIPTION
C  -------- --- ---- --- -----------
C
C  TSEC50    1  R*8   I  TIME ASSOCIATED WITH THE MEAN OF DATE
C                        COORDINATE SYSTEM. I.E., ROTATION IS BETWEEN
C                        THE MEAN OF 1950.0 AND THE MEAN OF TSEC50
C                        COORDINATE SYSTEMS. IN SECONDS SINCE
C                        0.0 HOURS, 1/1/50, ET. ACCURACY IS NOT HIGH
C                        ENOUGH TO WARRANT ET-UT CORRECTION. TIME MAY BE
C                        UT IF IT IS MORE CONVENIENT.
C
C  KFLAG     1  I*4   I  FLAG TELLING ROUTINE DIRECTION OF ROTATION.
C                        = NEGATIVE, FROM MEAN OF DATE TO MEAN OF 1950.0
C                        = OTHERWISE, FROM MEAN OF 1950.0 TO MEAN OF DATE.
C
C  ROTMTX   3,3 R*8   O  THE ROTATION MATRIX. ROTATION IS DONE
C                        EXTERNALLY BY:
C                         VEC2(I)=ROTMTX(I,J)*VEC1(J), SUMMING ON J.
C                         IF KFLAG IS NEGATIVE, VEC1 IS MEAN OF DATE.
C                         OTHERWISE, VEC1 IS MEAN OF 1950.0
C
C**********************************************************************
C
C  BY C PETRUZZO, 12/82.
C   MODIFIED.........
C           7/83. CJP. CHANGED DAYCENT VALUE.
C
C**********************************************************************
C
C  >> NOTE: FORMULAE AND CONSTANTS WERE TAKEN FROM THE EXPLANATORY
C           SUPPLEMENT TO THE ASTRONOMICAL EPHEMERIS AND THE AMERICAN
C           EPHEMERIS AND NAUTICAL ALMANAC, 1974(=1961 + AMMENDMENTS).
C           IN THE GSFC LIBRARY, NUMBER= QB8.U5, E96.
C
      REAL*8  A(3,3),ROTMTX(3,3),SECDAY/86400.D0/,TLAST/-1.D30/
C     SET NUMBER OF DAYS IN A TROPICAL CENTURY. SUPPLEMENT, PAGE 22.
      REAL*8  DAYCENT/36524.21988D0/
C
      IF(TSEC50.EQ.TLAST) GO TO 800
C
      TLAST=TSEC50
      T=TSEC50/SECDAY+0.077D0     ! DAYS SINCE BESSELIAN YEAR 1950.0
      T1 = T/DAYCENT              ! TROPICAL CENTURIES SINCE BYR 1950.0
      T2 = T1 * T1
      T3 = T2 * T1
C
C     GET ROTATION MATRIX FOR MEAN OF 1950.0 TO MEAN OF DATE. PAGE 34.
      A(1,1) = 1.D0                   - 29696.D-8 * T2  -  13.D-8 * T3
      A(1,2) =      -2234941.D-8 * T1 -   676.D-8 * T2  + 221.D-8 * T3
      A(1,3) =      - 971690.D-8 * T1 +   207.D-8 * T2  +  96.D-8 * T3
      A(2,1) = -A(1,2)
      A(2,2) = 1.D0                   - 24975.D-8 * T2  -  15.D-8 * T3
      A(2,3) =                        - 10858.D-8 * T2
      A(3,1) = -A(1,3)
      A(3,2) =  A(2,3)
      A(3,3) = 1.D0                   -  4721.D-8 * T2
C
  800 CONTINUE
      IF(KFLAG.LT.0) THEN
C       MEAN OF DATE TO MEAN OF 1950.
        DO 20 I=1,3
        DO 20 J=1,3
   20   ROTMTX(I,J)=A(J,I)
      ELSE
C       MEAN OF 1950.0 TO MEAN OF DATE.
        DO 30 I=1,3
        DO 30 J=1,3
   30   ROTMTX(I,J)=A(I,J)
        END IF
C
      RETURN
      END

ext: for


module zgrid

  implicit none

  public :: init_zgrid, finish_zgrid
  public :: nzed, nzgrid, nperiod, ntubes
  public :: nztot, nz2pi
  public :: zed
  public :: delzed
  public :: zed_equal_arc
  public :: get_total_arc_length
  public :: get_arc_length_grid
  public :: shat_zero
  public :: boundary_option_switch
  public :: boundary_option_zero
  public :: boundary_option_self_periodic
  public :: boundary_option_linked

  private

  integer :: nzed, nzgrid, nztot, nz2pi
  integer :: nperiod, ntubes
  logical :: zed_equal_arc
  real :: shat_zero
  real, dimension (:), allocatable :: zed, delzed

  integer :: boundary_option_switch
  integer, parameter :: boundary_option_zero = 1, &
       boundary_option_self_periodic = 2, &
       boundary_option_linked = 3

  logical :: zgridinit = .false.

contains

  subroutine init_zgrid

    use mp, only: proc0
    use constants, only: pi

    implicit none

    integer :: i

    if (zgridinit) return
    zgridinit = .true.

    if (proc0) then
       call read_parameters
    end if
    call broadcast_parameters

    if (.not.allocated(zed)) allocate (zed(-nzgrid:nzgrid))
    if (.not.allocated(delzed)) allocate (delzed(-nzgrid:nzgrid))

    zed = (/ (i*pi/real(nzed/2), i=-nzgrid, nzgrid ) /)
    delzed(:nzgrid-1) = zed(-nzgrid+1:) - zed(:nzgrid-1)
    delzed(nzgrid) = delzed(-nzgrid)

    nztot = 2*nzgrid+1
    ! number of zed in a 2*pi segment, including points at +/- pi
    nz2pi = 2*(nzed/2)+1

  end subroutine init_zgrid

  subroutine read_parameters

    use file_utils, only: input_unit_exist, error_unit
    use text_options, only: text_option, get_option_value
    use physics_flags, only: full_flux_surface

    implicit none

    integer :: in_file, ierr
    logical :: exist

    type (text_option), dimension (6), parameter :: boundaryopts = &
         (/ text_option('default', boundary_option_zero), &
            text_option('zero', boundary_option_zero), &
            text_option('unconnected', boundary_option_zero), &
            text_option('self-periodic', boundary_option_self_periodic), &
            text_option('periodic', boundary_option_self_periodic), &
            text_option('linked', boundary_option_linked) /)
    character(20) :: boundary_option

    namelist /zgrid_parameters/ nzed, nperiod, ntubes, &
         shat_zero, boundary_option, zed_equal_arc

    nzed = 24
    nperiod = 1
    ntubes = 1
    boundary_option = 'default'
    ! if zed_equal_arc = T, then zed is chosen to be arc length
    ! if zed_equal_arc = F, then zed is poloidal (axisymmetric)
    ! or zeta (toroidal) angle
    zed_equal_arc = .false.
    ! set minimum shat value below which we assume
    ! periodic BC
    shat_zero = 1.e-5

    in_file = input_unit_exist("zgrid_parameters", exist)
    if (exist) read (unit=in_file, nml=zgrid_parameters)

    ierr = error_unit()
    call get_option_value &
         (boundary_option, boundaryopts, boundary_option_switch, &
         ierr, "boundary_option in dist_fn_knobs")

    ! note that boundary_option may be changed to self-periodic later
    ! if magnetic shear is smaller than shat_zero
    ! cannot do this here as magnetic shear has yet to be input

    nzgrid = nzed/2 + (nperiod-1)*nzed

    ! force use of equal arc grid to ensure gradpar alpha-independent
    ! necessary to obtain efficient numerical solution of parallel streaming
    if (full_flux_surface) zed_equal_arc = .true.

  end subroutine read_parameters

  subroutine broadcast_parameters

    use mp, only: broadcast

    implicit none

    call broadcast (nzed)
    call broadcast (nzgrid)
    call broadcast (nperiod)
    call broadcast (ntubes)
    call broadcast (zed_equal_arc)
    call broadcast (shat_zero)
    call broadcast (boundary_option_switch)

  end subroutine broadcast_parameters

  subroutine finish_zgrid

    implicit none

    if (allocated(zed)) deallocate (zed)
    if (allocated(delzed)) deallocate (delzed)

    zgridinit = .false.

  end subroutine finish_zgrid

  subroutine get_total_arc_length (nz, gp, dz, length)

    implicit none

    integer, intent (in) :: nz
    real, dimension (-nz:), intent (in) :: gp
    real, intent (in) :: dz
    real, intent (out) :: length

    call integrate_zed (nz, dz, 1./gp, length)

  end subroutine get_total_arc_length

  subroutine get_arc_length_grid (nz_max, nzext_max, zboundary, gp, dz, zarc)

    implicit none

    integer, intent (in) :: nz_max, nzext_max
    real, intent (in) :: zboundary
    real, dimension (-nzext_max:), intent (in) :: gp
    real, intent (in) :: dz
    real, dimension (-nzext_max:), intent (out) :: zarc

    integer :: iz

    zarc(-nz_max) = zboundary
    if (nz_max /= nzext_max) then
       do iz = -nzext_max, -nz_max-1
          call integrate_zed (nzext_max, dz, 1./gp(iz:-nz_max), zarc(iz))
          zarc(iz) = zarc(-nz_max) - zarc(iz)
       end do
    end if
    do iz = -nz_max+1, nzext_max
       call integrate_zed (nz_max, dz, 1./gp(-nz_max:iz), zarc(iz))
       zarc(iz) = zarc(-nz_max) + zarc(iz)
    end do

  end subroutine get_arc_length_grid

  ! trapezoidal rule to integrate in zed
  subroutine integrate_zed (nz, dz, f, intf)

    implicit none

    integer, intent (in) :: nz
    real, intent (in) :: dz
    real, dimension (-nz:), intent (in) :: f
    real, intent (out) :: intf

    integer :: iz, iz_max

    iz_max = -nz + size(f) - 1
    intf = 0.
    do iz = -nz+1, iz_max
       intf = intf + dz*(f(iz-1)+f(iz))
    end do
    intf = 0.5*intf

  end subroutine integrate_zed

end module zgrid


text:
This is a toy, book, and game store focusing on natural and Davis Waldorf School Waldorfinspired products encouraging childrens natural creativity and imagination.

As of August 2007, this store has closed its doors. No note on the door sheds any light as to its closing, but the room is empty of all inventory. The space has been the site of at least three games/toy stores in the past: The Game Preserve (Derrick Bang), The End Zone (run by Vincentor Ian?) where kids/adults could buy and play card games like Magic, and now Aesops. Sad.

Recently I went in to buy a gift for my nephew and ended up with an awesome set of Anker architecture blocks as well as some interesting logic puzzles. I noticed a cool collection of math and science toys happening.  DavidPonce

20070112 22:49:01 nbsp I love the products but cant afford most of them. A kite cost me $25! It was a nice kite though. Nothing like the cheap stuff I grew up with. It better last a long time... RoseNylund

20070806 12:44:40 nbsp I was sad to see Aesops closed last week...is this now a departed business? Users/ScottLay

20070821 21:36:45 nbsp Aesops Room was having hard times due to a slow economy, and the owners couldnt keep it going as a storefront business.  It is still open as an internetbased toy store, http://www.aesopsroom.com .  Sad to see it go.  Rememberif you value your downtown shops, make sure you buy stuff from them!   Users/TedBuehler

20090826 21:36:45 nbsp I was the original creator and owner of Aesops Room. Having began, grown and operated it for about seven years at that point, I was ready to move on. After approximately a year under new ownership, it was disheartening to see Aesops close. This was the second business I have sold, the first still alive and thriving. I now dedicate myself to helping others grow their business and meet their goals through http://lifecompassnow.com Coaching and http://carmenisais.com Business Consultation.  Users/Cisais

20090828 21:39:49 nbsp I only worked here a short time before they closed, but loved it. Carmen was great. Very sad to see the new owners couldnt hack it.  Users/EricaMacGregor

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text: #redirect Fiddlers Green


text:
Describe Users/erindunning here.

20110923 13:41:31 nbsp Hi, Erin, and welcome to the wiki!  Thanks for your updates to the Bikram Yoga Davis page.  You wondered who was in charge of the Davis wiki  the answer is, we all are!  That is, all of the members of the Davis community who care to participate.  So, yes, when you see something wrong, please update it.  You might also want to read the Welcome to the wiki/Business Owner page.  It tries to clarify the roles of business owners and and the community in editing pages about Davis businesses.  I hope you stick around to make other edits! Users/CovertProfessor

20110923 14:05:48 nbsp Hi Erin Users/StevenDaubert

20110923 14:21:42 nbsp I hacked the Legend of Yona page. Users/jefftolentino
    Watch out, now. In case she unsheathes her hacker/peacemaker as she looks prone to do on general principles, anyways. Users/RaoulDuke

20111005 11:03:57 nbsp hey mr. covert... i used hack loosely to be funny. i figure as the studio owner i had the right to put in accurate information!  Users/erindunning
  Absolutely!  cp

20111005 11:05:33 nbsp jeff!  theres a yona page??  Users/erindunning
 The Saga of Yona Of course!  Users/WilliamLewis

20111005 13:29:47 nbsp Hmm, I am not sure if I can take the heat  literally!  Im still adjusting to the Davis climate as it is. Users/CovertProfessor


text:
module rp1_advection
  use iso_c_binding
  implicit none

  type :: advection_context
     real(c_double) :: u
  end type advection_context

contains

  subroutine new(u, context) bind(C, name='rp1_advection_new')
    real(c_double), intent(in)       :: u
    type(c_ptr), intent(out)         :: context

    type(advection_context), pointer :: param

    allocate(param)
    param%u = u
    context = c_loc(param)
  end subroutine new

  subroutine delete(context) bind(C, name='rp1_advection_delete')
    type(c_ptr), intent(in)          :: context

    type(advection_context), pointer :: param

    call c_f_pointer(context, param)
    deallocate(param)
  end subroutine delete

  ! only 1 wave, no need for an mwave parameter
  ! only 1 equation, no need for an meqn parameter

  subroutine advection(mx, mbc, ql, qr, s, wave, amdq, apdq, context) &
       bind(C, name='rp1_advection_c')
    ! C interfacable version of Riemann solver for advection
    !
    !
    ! inputs
    ! * ql  - data at left edge of each cell
    ! * qr  - data at right edge of each cell
    ! * mx  - number of physical grid points
    ! * mbc - number of boundary layer grid points
    ! outputs
    ! * s
    ! * wave
    ! * amdq
    ! * apdq
    ! context
    ! * u - constant advection speed
    !

    ! Solve Riemann problems for the 1D advection equation q_t + u*q_x = 0.
    ! for constant advection velocity u, passed in common block.
    !
    ! The advection speed u is passed in the derived type solver data
    ! On input, ql contains the state vector at the left edge of each cell
    !           qr contains the state vector at the right edge of each cell
    ! On output, wave contains the waves,
    !            s the speeds,
    !            amdq the  left-going flux difference  A^- \Delta q
    !            apdq the right-going flux difference  A^+ \Delta q
    !
    ! Note that the i'th Riemann problem has left state qr(i-1,:) and
    !                                    right state ql(i,:) From the
    !                                    basic clawpack routine step1, rp
    !                                    is called with ql = qr = q.

    integer(c_int), intent(in), value :: mx, mbc
    real(c_double), intent(in) :: &
         ql(1, 1-mbc:mx+mbc), &
         qr(1, 1-mbc:mx+mbc)
    type(c_ptr), intent(in), value :: context

    real(c_double), intent(out) :: &
         s(1, 1-mbc:mx+mbc), &
         wave(1, 1, 1-mbc:mx+mbc), &
         apdq(1, 1-mbc:mx+mbc), &
         amdq(1, 1-mbc:mx+mbc)

    type(advection_context), pointer :: param
    real(c_double) :: u
    integer :: i

    call c_f_pointer(context, param)
    u = param%u

    do i=2-mbc,mx+mbc
       ! Compute the wave
       wave(1,1,i) = ql(1,i) - qr(1,i-1)

       ! Speeds are all the same
       s(1,i) = u

       ! Set fluctuations
       amdq(1,i) = dmin1(u, 0.d0) * wave(1,1,i)
       apdq(1,i) = dmax1(u, 0.d0) * wave(1,1,i)
    end do

  end subroutine advection

end module rp1_advection

max_stars_repo_name: gyanz/pyclaw

max_issues_repo_name: gyanz/pyclaw

max_forks_repo_name: gyanz/pyclaw


text:
C **********************************************************************
C                                                  Version 2.0
C     Generator of e+e- pairs produced in            1.10.2002
C                  PbPb collisions at LHC
C
C     Copyright (c) 2002
C     Authors:          Kai Hencken	     <k.hencken@unibas.ch>
C                       Yuri Kharlov	     <Yuri.Kharlov@cern.ch>
C                       Serguei Sadovsky     <sadovsky@mx.ihep.su>,
C                                            <Serguei.Sadovski@cern.ch>
C
C Permission to use, copy and distribute this software and its
C documentation strictly for non-commercial purposes is hereby granted
C without fee, provided that the above copyright notice appears in all
C copies and that both the copyright notice and this permission notice
C appear in the supporting documentation. The authors make no claims
C about the suitability of this software for any purpose. It is
C provided "as is" without express or implied warranty.
C Any change of the code should be submitted to the authors
C
C     Long write up:   K.Hencken,Yu.Kharlov,S.Sadovsky
C                      Internal ALICE Note 2002-27
C
C **********************************************************************
      subroutine ee_init(Ymin,Ymax,PTmin,PTmax)
C-----------------------------------------------------------------------
C    Generator initialisation
C
C    Input variables:  Ymin      - minimal value of rapidity  \
C                      Ymax      - maximal value of rapidity   | of kinematics
C                      PTmin     - Pt minimum in MeV/c	       | range
C                      PTmax     - Pt maximum in MeV/c        /
C-----------------------------------------------------------------------
      implicit real*8 (A-H,O-Z)
      external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi,DsdXX,DsdYY

      common/parPh / parPhi(7)
      common/parYp / parYpY(6)
      common/parYm / parYmY(6)
      common/parXp / parXpX(5)
      common/parXm / parXmX(4)

      common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt
      common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax,
     &               Ymed1,Ymed2,sgmY1,sgmY2,XYsect,
     &               Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3,
     &               XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed,
     &               Sgm1,AnorX,Icase

      real*8 mass

      data pi     / 3.141 592 653 589 793 238 462 643d00 /

C
C  Exact differential Cross section:
C - - - - - - - - - - - - - - - - - - -
C
      gm = 2750.0d0             ! Pb gamma factor at the LHC
      mass = 0.5109991d0        ! electron mass
      call initdiffcross (gm,mass)
C
C  Cross sections:
C - - - - - - - - - -
      Nevnt  = 0
      Xsect2 = 0.
      Dsect2 = 0.

      Nsd = 1
      Npt = 25                  !  7,25,64
      Eps = 0.00005

      Xmin=    Dlog10(Ptmin)
      Xmax=    Dlog10(Ptmax)
      write(*,*) ' Kinematical limits:'
      write(*,*) '          Xmin,Xmax=',Xmin,Xmax
      write(*,*) '          Ymin,Ymax=',Ymin,Ymax
      write(*,*)

      Xsec1  = Dtrint(DsdXX,Nsd,Npt,Eps,Xmin,Xmin,Xmin,Xmax,Xmax,Xmin)
      Xsec2  = Dtrint(DsdXX,Nsd,Npt,Eps,Xmax,Xmax,Xmin,Xmax,Xmax,Xmin)
      XsecX  = Xsec1+Xsec2

      Ysec1  = Dtrint(DsdYY,Nsd,Npt,Eps,Ymin,Ymin,Ymin,Ymax,Ymax,Ymin)
      Ysec2  = Dtrint(DsdYY,Nsd,Npt,Eps,Ymax,Ymax,Ymin,Ymax,Ymax,Ymin)
      XsecY  = Ysec1+Ysec2

      IF (Xsec1*Xsec2.le.0.or.Ysec1*Ysec2.le.0.) then
         write(*,*) ' Error: insufficient accuracy of XY-cross sections'
         write(*,*) ' Xsec1,Xsec2,Xsec=', Xsec1,Xsec2,XsecX
         write(*,*) ' Ysec1,Ysec2,Ysec=', Ysec1,Ysec2,XsecY
      endif

      XYsect = XsecX*XsecY
      write(*,*) ' Xsections: Xsec1,Xsec2,XsecX=', Xsec1,Xsec2,XsecX
      write(*,*) '            Ysec1,Ysec2,XsecY=', Ysec1,Ysec2,XsecY
C-    write(*,*) '                  XsecX*YsecY=', XYsect
C
      XYsect = 2371.5239*(82./137.035)**4*XYsect ! Normalization factor
      write(*,*) ' Normalized:      XsecX*YsecY=', XYsect
      write(*,*)
C
C-		   XsecPhi  = Dgauss(DsdPhi,0.0d0,  pi, eps)
C
C  Rapidity initialisation:
C - - - - - - - - - - - - - -
      if (Ymin.ge.Ymax) then
         write(*,*) 'Wrong values of Ymin,Ymax:',Ymin,Ymax
         stop
      endif
CYpY-
C
      YpYmin = 2.*Ymin
      YpYmax = 2.*Ymax
      Yp0    = 0.
      if (YpYmin*YpYmax.gt.0.) then
         if(YpYmin.gt.0.) Yp0 = YpYmin
         if(YpYmin.le.0.) Yp0 = YpYmax
      endif
      WYpYmax = DsdYpY(Yp0)
C
C-    write(*,*)' YpY: YpYmin,YpYmax,WYpYmax =',YpYmin,YpYmax,WYpYmax

CYmY-
C
      YmYmin = Ymin-Ymax
      YmYmax = Ymax-Ymin
C-    write(*,*) ' YmY:  YmYmin,YmYmax=',YmYmin,YmYmax
C
      Ymed1  = 0.18
      Ymed2  = 4.00

      sgmY1  = 1./dsqrt(2.*parYmY(2))
      sgmY2  = 1./dsqrt(2.*parYmY(4))
C
      eps = 0.000001
      Gaus1 = Dgauss(DsdYmY ,0.0d0, Ymed1 , eps)
      Gaus2 = Dgauss(DsdYmY ,Ymed1, Ymed2 , eps)
      Expon = Dgauss(DsdYmY ,Ymed2, YmYmax, eps)
      Summ  = Gaus1+Gaus2+Expon
      Gaus2 =(Gaus1+Gaus2)/Summ
      Gaus1 = Gaus1/Summ
      Expon = 1.
C
C     Phi initialisation:
C   - - - - - - - - - - -
      Exp0 = parPhi(1)*pi
      Exp1 = parPhi(2)*(1.-dexp(-parPhi(3)*pi))/parPhi(3)
      Exp2 = parPhi(4)*(1.-dexp(-parPhi(5)*pi))/parPhi(5)
      Exp3 = parPhi(6)*(1.-dexp(-parPhi(7)*pi))/parPhi(7)
      Summ = Exp0+Exp1+Exp2+Exp3
C
      Exp0 = (Exp0+Exp1+Exp2+Exp3)/Summ
      Exp1 = (Exp1+Exp2+Exp3)/Summ
      Exp2 = (Exp2+Exp3)/Summ
      Exp3 =  Exp3/Summ
C-    write(*,*) 'Exp0,Exp1,Exp2,Exp3=',Exp0,Exp1,Exp2,Exp3
C
C  Pt initialisation:
C - - - - - - - - - - -
C
      XmXmin = Xmin-Xmax
      XmXmax = Xmax-Xmin
C
      XpXmin = 2.*Xmin
      XpXmax = 2.*Xmax
C
CXmX-
      Exmx1 = parXmX(1)*(1.-dexp(-parXmX(2)*dabs(XmXmax)))/parXmX(2)
      Exmx2 = parXmX(3)*(1.-dexp(-parXmX(4)*dabs(XmXmax)))/parXmX(4)
      Exmx1 = Exmx1/(Exmx1+Exmx2)
      Exmx2 = 1.
C
C-    write(*,*) ' XpX:  XpXmin,XpXmax=',XpXmin,XpXmax
C-    write(*,*) ' XmX:  XmXmin,XmXmax=',XmXmin,XmXmax
C-    write(*,*) ' XmX:  Exmx1 ,Exmx2 =',Exmx1 ,Exmx2
C
CXpX-
      Icase = 2                 !   Gausses and Exponent
      Xmed  = 0.6
      if (XpXmax.lt.Xmed) then
         Icase = 1              !   Gauss only
         Xmed  = XpXmax
      endif
C
      if (XpXmin.gt.Xmed) then
         Icase = 3              !   Exponent only
         Xmed  = XpXmin
      endif
C
      if (Icase.eq.2)     then  !   Gauss and Exponent
         eps = 0.000001
         Gauss = Dgauss(DsdXpX ,XpXmin, Xmed  , eps)
         Expon = Dgauss(DsdXpX ,Xmed  , XpXmax, eps)
         Summ  = Gauss+Expon
         Gauss = Gauss/Summ
         Expon = Expon/Summ
C-       write(*,*) 'XpX, Case 2: Gauss,Expon=',Gauss,Expon
      endif                     ! Icase = 2
C
C
      if (Icase.eq.1.or.Icase.eq.2) then !   Gausses only and Icase=2

         Sgm1 =       1./dsqrt(2.*parXpX(2))

C-       write(*,*) 'XpX, Case 1:  Sgm1=',Sgm1
      endif
C
      if (Icase.eq.3.or.Icase.eq.2) then !    Exponent only and Icase=2
C
         AnorX = 1.-dexp(-parXpX(5)*(XpXmax-Xmed))
C
C-       write(*,*) 'XpX, Case 3:  AnorX=',AnorX
      endif

      return
      end
*
*=======================================================================
      subroutine ee_event(Ymin,Ymax,PTmin,PTmax,
     +                    Ye,Yp,Xe,Xp,Phi,Wtm2)
C------------------------------------------------------------------------------
C    Produce one event
C
C    Input variables:  Ymin     - minimal value of rapidity  \
C                      Ymax     - maximal value of rapidity   | of kinematics
C                      PTmin    - Pt minimum in MeV/c	      | range
C                      PTmax    - Pt maximum in MeV/c	     /
C    Output variables: Ye,Yp    - rapidity of produced e- or e+
C                      Xe,Xp    - log10(pt) of produced e- or e+, pt in MeV/c
C                      Phi      - azymuth angle between e- and e+
C                      Wtm2     - event weight. The sum of these event
C                                 weights, divided by the total number
C                                 of generated events, gives the integral
C                                 cross section of the process of e+e- pair
C  		                  production in the above mentioned kinematics
C                                 range.
C		                  Sum of the selected event weights, divided
C                                 by the total number of generated events,
C                                 gives the integral cross section corresponded
C                                 to the set of selected events
C------------------------------------------------------------------------------

      implicit real*8 (A-H,O-Z)
      common/parPh / parPhi(7)
      common/parYp / parYpY(6)
      common/parYm / parYmY(6)
      common/parXp / parXpX(5)
      common/parXm / parXmX(4)
      common/eevent/ Xsect2,Dsect2, Xsecttot,Dsecttot, Nevnt
      common/eepars/ Xmin,Xmax,YpYmin,YpYmax,WYpYmax,YmYmin,YmYmax,
     &               Ymed1,Ymed2,sgmY1,sgmY2,XYsect,
     &               Gaus1,Gaus2,Gauss,Exp1,Exp2,Exp3,
     &               XmXmin,XmXmax,XpXmin,XpXmax,Exmx1,Exmx2,Xmed,
     &               Sgm1,AnorX,Icase
      external DsdYpY,DsdYmY,DsdXpX,DsdXmX,DsdPhi
      data pi     / 3.141 592 653 589 793 238 462 643d00 /

C  Rapidity distributions:
C
 10   YpY = YpYmin + (YpYmax-YpYmin)*eernd(0)
      Wt  = DsdYpY(YpY)
      if (Wt.lt.eernd(0)*WYpYmax) go to 10
C
CYmY-
      r1 = eernd(0)
      if  (r1.lt.Gaus1) then
 13      call rnormlEP(rn1)
         YmY = sgmY1*rn1
         if (dabs(YmY).gt.Ymed1) go to 13
      else if  (r1.lt.Gaus2) then
 15      call rnormlEP(rn2)
         YmY = sgmY2*rn2
         if (dabs(YmY).lt.Ymed1) go to 15
         if (dabs(YmY).gt.Ymed2) go to 15
      else
 20      r2 = eernd(0)
         YmY=-Dlog(r2)/parYmY(6) + Ymed2
         if (YmY.gt.YmYmax) go to 20
         r2 = eernd(0)
         if (r2.lt.0.5) YmY =-YmY
      endif
C
      Ye = 0.5*(YpY+YmY)
      Yp = 0.5*(YpY-YmY)

      if (Ye.lt.Ymin.or.Ye.gt.Ymax) go to 10
      if (Yp.lt.Ymin.or.Yp.gt.Ymax) go to 10

C
C Azimuthal angle:
C
 30   r1 = eernd(0)
      if (r1.lt.Exp3) then
 31      r2 = eernd(0)
         Phi=-Dlog(r2)/parPhi(7)
         if (Phi.gt.pi) go to 31
         go to 50
      else if (r1.lt.Exp2) then
 32      r2 = eernd(0)
         Phi=-Dlog(r2)/parPhi(5)
         if (Phi.gt.pi) go to 32
         go to 50
      else if (r1.lt.Exp1) then
 33      r2 = eernd(0)
         Phi=-Dlog(r2)/parPhi(3)
         if (Phi.gt.pi) go to 33
         go to 50
      else
         Phi = pi*eernd(0)
      endif
C
 50   if (eernd(0).gt.0.5) Phi =-Phi
      Phi = Phi+pi
      if (Phi.lt.0.03.or.Phi.gt.2.*pi-0.03) go to 30
C
C Transverse momentums:
C
CXpX-

  60  Jcase = Icase             !       Jcase = Icase, if Icase.ne.2
C
      if (Icase.eq.2) then      !       Gausses and Exponent
         Jcase = 3
         if (eernd(0).lt.Gauss) Jcase = 1
      endif                     ! of Icase 2
C
      if (Jcase.eq.1) then      !       Gauss only
 65      call rnormlEP(rn1)
         XpX=Sgm1*rn1-parXpX(3)
         if (XpX.lt.XpXmin.or.XpX.gt.Xmed) go to 65
      endif                     ! of Jcase 1
C
      if (Jcase.eq.3) then      !       Exponent only
 70      r1  = eernd(0)
         XpX =-Dlog( 1.- r1*AnorX)/parXpX(5)+Xmed
      endif                     ! of Jcase 3

CXmX-
      r1 = eernd(0)
      if (r1.lt.Exmx1) then
 81      r2 = eernd(0)
         XmX=-Dlog(r2)/parXmX(2)
         if (XmX.gt.dabs(XmXmax)) go to 81
      else
 82      r2 = eernd(0)
         XmX=-Dlog(r2)/parXmX(4)
         if (XmX.gt.dabs(XmXmax)) go to 82
      endif
C
 85   if (eernd(0).gt.0.5) XmX =-XmX
C
      Xe = 0.5*(XpX+XmX)
      Xp = 0.5*(XpX-XmX)

      if (Xe.lt.Xmin.or.Xe.gt.Xmax) go to 60
      if (Xp.lt.Xmin.or.Xp.gt.Xmax) go to 60
C
C --- Good event:
C
      Pte = 10.d0**Xe
      Ptp = 10.d0**Xp
C
C Event weight if the events would generate uniformly
C
      Wt = DsdYpY(YpY)*DsdYmY(YmY)*DsdXpX(Xpx)*DsdXmX(XmX)*DsdPhi(Phi)
      Wt = 2.2706950/Wt
C
      Wtm2=(Pte*Ptp)*Wt          ! Transformation factor
C
C --- Exact diff. cross section (Adrian Alscher 1997, Kai Hencken, May '98):
C
      call Diffcross(Ptp,Yp,Pte,Ye,Phi,Dsigma)
      Wtm2= XYsect*Dsigma*(Pte*Ptp)*Wtm2
C
      Nevnt  = Nevnt  + 1
      Xsect2 = Xsect2 + Wtm2
      Dsect2 = Dsect2 + Wtm2*Wtm2
C
C     Calculate cross section and its error accumulated so far
C
      Xsecttot = Xsect2/Nevnt
      Dsecttot = dsqrt(Dsect2)/Nevnt
C
      return
      end
c
C=================================================================
      Double precision function DsdYpY(Y)
c
c     Y = Yp+Ye
c - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
      common/parYp/   parYpY(6)
C-    data parYpY /  -48.584, 0.11403E-01, 40.649, 0.38920E-04,
C-   +                                     40.283, 0.19238E-01 /
      data parYpY /  -4.8584, 0.11403E-01, 4.0649, 0.38920E-04,
     +                                     4.0283, 0.19238E-01 /
c
      DsdYpY = parYpY(1)*dexp(-parYpY(2)*Y**2)
     +       + parYpY(3)*dexp(-parYpY(4)*Y**4)
     +       + parYpY(5)*dexp(-parYpY(6)*Y**2)
      return
      end
C
      Double precision function DsdYmY(Y)
c
c     Y = Yp-Ye
c - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
      common/parYm/   parYmY(6)
C-    data parYmY /   180.,  8., 174.38, 0.17867, 2418.8, 1.3097 /
      data parYmY /   1.80,  8., 1.7438, 0.17867, 24.188, 1.3097 /
C
      if (abs(Y).lt.0.18) then
C1-      if (abs(Y).lt.0.00) then
         DsdYmY  = parYmY(1)*dexp(-parYmY(2)*abs(Y)**2)
      else if (abs(Y).lt.4.00) then
         DsdYmY  = parYmY(3)*dexp(-parYmY(4)*abs(Y)**2)
      else
         DsdYmY  = parYmY(5)*dexp(-parYmY(6)*abs(Y))
      endif

      return
      end
C
      Double precision function DsdYY(Yp,Ye)
c - - - - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
c
      YpY   = Yp+Ye
      YmY   = Yp-Ye
      DsdYY = DsdYpY(YpY)*DsdYmY(YmY)
      return
      end
C
C====================================================================
C
      Double precision function DsdXpX(X)
c
c     X = Xp+Xe
c - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
      common/parXp/   parXpX(5)

C-    data parXpX /   83.668, 1.2004, 0.47225, 96.951, 2.3814 /
      data parXpX /   8.3668, 1.2004, 0.47225, 9.6951, 2.3814 /
c
      if (X.lt.0.6)then
         DsdXpX = parXpX(1)*dexp(-parXpX(2)*(X+parXpX(3))**2)
      else
         DsdXpX = parXpX(4)*dexp(-parXpX(5)*X)
      endif
      return
      end
C
      Double precision function DsdXmX(X)
c
c     X = Xp-Xe
c - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
      common/parXm/   parXmX(4)
c-    data parXmX /   950.38, 39.040, 194.92, 3.5660 /
      data parXmX /   9.5038, 39.040, 1.9492, 3.5660 /
c

      DsdXmX = parXmX(1)*dexp(-parXmX(2)*abs(X))
     +       + parXmX(3)*dexp(-parXmX(4)*abs(X))

      return
      end
C
      Double precision function DsdXX(Xp,Xe)
c - - - - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
c
      XpX   = Xp+Xe
      XmX   = Xp-Xe
      DsdXX = DsdXpX(XpX)*DsdXmX(XmX)
      return
      end
C====================================================================
C
      Double precision function DsdPhi(Phi)  ! Phi distribution
c
c     Phi - athimutal angle between e+ and e- (in radians)
c - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit real*8 (A-H,O-Z)
      common/parPh/ parPhi(7)
C                                                         !  For 10**7 events
      data parPhi / 5.4825, 226.00, 11.602,      68.173,  !  accuracy= 0.87 %
     +			    1.9509, 0.11864E+04, 109.61 / !  N(Wt>20)= 573

C
      data pi     / 3.141 592 653 589 793 238 462 643d00/
C
      DsdPhi  = parPhi(1)
     +        + parPhi(2)*dexp(-parPhi(3)*dabs(Phi-pi))
     +        + parPhi(4)*dexp(-parPhi(5)*dabs(Phi-pi))
     +        + parPhi(6)*dexp(-parPhi(7)*dabs(Phi-pi))
      return
      end
C====================================================================
      SUBROUTINE rnormlEP(rnd)
*     Random generator of normal distribution
      implicit real*8 (A-H,O-Z)
      PARAMETER (pi=3.141 592 653 589 793 238 462 643d00)
      PARAMETER (pi2=2.*pi)
C
      u1 = eernd(0)
      u2 = eernd(0)
      p1 = pi2*u1
      p2 = dsqrt(-2.*dlog(u2))
      rnd= dcos(p1)*p2
      RETURN
      END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot