      real function simplem(x) 


      vector pars(9) 

      vector den(1)  

      vector opar(10) 

      real x 
      
      real th1, th2, th3
      real s1, s2, s3 
      real c1, c2, c3
     
      real delta
      real sd, cd
      real jcp 
      real icp 
      real dm12 
      real dm23 
      real dm31 
      real ph1, ph2, ph3 
      real ord
      real anu

      real term1, term2, term3, term4 
      real termcp 

      real mtrm1, mtrm2, mtrm3, mtrm4 

      real density, ahat 

      ord = pars(8) 
      anu = pars(9) 

      dm12 = pars(1)
      dm31 = ord*pars(2)
      dm32 = ord*pars(2)-pars(1)  
      th1 = pars(3) 
      th2 = pars(4)
      th3 = pars(5) 
      delta = anu*pars(6) 
      dist = pars(7) 


c      write(6,*) pars
      
      ph1 = 1.27*dm12*dist 
      ph2 = 1.27*dm32*dist*ord
      ph3 = 1.27*dm31*dist*ord  

      s1 = sin(th1)
      s2 = sin(th2) 
      s3 = sin(th3) 
      c1 = cos(th1)
      c2 = cos(th2)
      c3 = cos(th3) 
      sd = sin(delta)
      cd = cos(delta) 

      jcp = s1*s2*s3*c1*c2*c3*c3*sd
      icp = s1*s2*s3*c1*c2*c3*c3*cd


      density = den(1) 

      ahat  = anu*7.63e-5*density*x/dm31 


      opar(1) = ahat
      opar(2) = s2*s2*4.*s3*s3*c3*c3
      opar(3) = jcp 
      opar(4)  = icp 
      opar(5)  = (dm12/dm31)**2*c2*c2*4.*s1*s1*c1*c1
      opar(6) = dm12/dm31*8.*jcp


      mtrm1 = s2*s2*4.*s3*s3*c3*c3*(sin((ahat-1)*ph3/x))**2/
     +        (ahat-1)**2  

      mtrm2 = dm12/dm31*8.*jcp*sin(ph3/x)*sin(ahat*ph3/x)
      mtrm2 = mtrm2*sin((1-ahat)*ph3/x)/ahat/(1-ahat) 

      mtrm3 = dm12/dm31*8.*icp*cos(ph3/x)*sin(ahat*ph3/x)
      mtrm3 = mtrm3*sin((1-ahat)*ph3/x)/ahat/(1-ahat) 

      mtrm4 = (dm12/dm31)**2*c2*c2*4.*s1*s1*c1*c1*(sin(ahat*ph3/x))**2
      mtrm4 = mtrm4/ahat/ahat 

      simplem = mtrm1 + mtrm2 + mtrm3 + mtrm4                       

      return 
      end 

