hump.f

Go to the documentation of this file.

      subroutine hump(const,t,tp,nm,lp,rmu,thd,ixy,nmu)                       
c                                                                       
c***********************************************************************
      implicit real*8(a-h,o-z)                                          
      real*8 t(1801,4),am(2500),cosn(2500),sn(2500),p(4,4)               
      real*8 tp(1000,4),cthr(2500),thr(2500),dmu(2500),rmu(51),thd(1801) 
c***********************************************************************
c    program to compute pbar for each mup when mu=mup and delfi=0       
c                                                                       
      izia=0
c     write(6,*)'entering hump'
      pi=dacos(-1.0d0 )
      conv =pi/180.0d0                                      
c                                                                       
c      for debug purposes, const=alambd**2/(4*pi**2*qs)
      do 300 i=3,4                                                      
      do 300 j=1,2                                                      
      p(i,j)=0.0d0                                                      
  300 p(j,i)=0.0d0                                                      
      dp=0.10d0
      dt=0.10d0
      i1=lp/2                                                           
      fi1=i1                                                            
      ni1=fi1/dp
      ni1=ni1+1
      dp=fi1/dfloat(ni1)
      c=180-i1                                                          
      i2=nm/2                                                           
c     write(*,*)'i1,i2,ni1',i1,i2,ni1
      do 2 k=1,ni1
      fk=dfloat(k-1)*dp
      if(ixy)60,61,60                                                   
   60 phi=(c+fk+dp/2.0d0)*conv                                             
      go to 62                                                          
   61 phi=(fk+dp/2.0d0)*conv                                               
62     sn(k) = dsin(phi)                                                
2      cosn(k) = dcos(phi)                                              
c     write(6,600) sn                                                   
 600  format(1x,'sn',1p10e12.4)                                         
      deldeg=dt
      do 7 kl=1,i2                                                      
      deg=rmu(kl+1) - rmu(kl)                                           
      nodeg=deg/deldeg
      deldeg=deg/dfloat(nodeg)                                          
      nodp1=nodeg+1                                                     
      xx=deg/2.0d0                                                      
      yy=xx-deldeg/2.0d0                                                
c                                                                       
      thetd=(rmu(kl)+rmu(kl+1))/2.0d0                                   
       amup = dcos(thetd*conv)                                          
c     if(thetd .lt. 3.1d0) amup=1.0d0                                   
c     write(6,9988)rmu(kl),rmu(kl+1),deg,thetd,nodeg                    
9988  format('rmuk1,rmuk2,deg,thetd,nodeg',1p4e11.3,i4)                 
      do 200 i=1,nodp1                                                  
      fi=dfloat(i-1)*deldeg                                             
  200 cthr(i)=(thetd-xx+fi)*conv                                        
      do 201 i=1,nodeg                                                  
      fi=dfloat(i-1)*deldeg                                             
      thr(i)=(thetd-yy+fi)*conv                                         
       am(i) = dcos(thr(i))                                             
  201 dmu(i)=dcos(cthr(i))-dcos(cthr(i+1))                              
c     write(6,602)(cthr(iz)/conv,iz=1,nodp1)                            
c     write(6,603)(thr(iz)/conv,iz=1,nodeg)                             
602   format('cthr',1p6e11.3)                                           
603   format('thr ',1p6e11.3)                                           
       ddmu =dcos(rmu(kl)*conv) -dcos(rmu(kl+1)*conv)                   
      amups=amup**2                                                     
      do 40 i=1,4                                                       
      k=(kl-1)*4+i                                                      
      do 40 j=1,4                                                       
   40 tp(k,j)=0.0d0                                                     
      do 1 i=1,nodeg                                                    
      if(ixy)500,501,500                                                
  500 amu=-am(i)                                                        
      go to 502                                                         
  501 amu=am(i)                                                         
  502 do 90 k=1,4                                                       
      do 90 l=1,4                                                       
   90 p(k,l)=0.0d0                                                      
      amusq=amu**2                                                      
      amumu=amu*amup                                                    
       anunu = dsqrt( (1.0d0 - amups)*(1.0d0 - amusq))                  
      do 150 j=1,ni1                                                    
      siin=sn(j)                                                        
      coosn=cosn(j)                                                     
      copsi=anunu+amumu*coosn                                           
      sifi=copsi*coosn                                                  
      sfisq=siin**2                                                     
      cfisq=coosn**2                                                    
      sisq=copsi**2                                                     
      a=amup*coosn                                                      
      b=amup*copsi                                                      
      c=amu*coosn                                                       
      d=amu*copsi                                                       
      e=sfisq*(amusq+amups)                                             
      g=amumu*sfisq                                                     
      h=sfisq*(amusq-amups)                                             
      x=amumu+anunu*coosn                                               
12     yz= dacos(x)/conv                                                
      yzz=10.0d0*yz
      m=idint(yzz)+1
      if(m.eq.0)m=1
c181  continue                                                          
c182  continue                                                          
      ya=yz-thd(m)                                                      
      dthd=thd(m)-thd(m+1)                                              
      call xntpln(yz,thd(m),thd(m+1),t(m,1),t(m+1,1),t1)
      call xntpln(yz,thd(m),thd(m+1),t(m,2),t(m+1,2),t2)
      call xntpln(yz,thd(m),thd(m+1),t(m,3),t(m+1,3),t3)
      call xntpln(yz,thd(m),thd(m+1),t(m,4),t(m+1,4),t4)
c     if(kl.le.1)then
c     write( 6,403)m,t(m,1),t(m+1,2),t1,yz,thd(m)                    
c     endif
c
      p(3,3)=p(3,3)+(sifi-g)*(t1+t2)+(cfisq+sisq-e)*t3                  
      p(2,2)=p(2,2)+sisq*t1+cfisq*t2+2.0d0 *sifi * t3                   
      p(1,2)=p(1,2)+sfisq*(amups*t1+amusq*t2+2.0d0 * amumu *t3 )        
      p(2,1)=p(2,1)+sfisq*(amups*t2+amusq*t1+2.0d0 * amumu *t3)         
      p(1,1)=p(1,1)+cfisq*t1+sisq*t2+2.0d0 * sifi *t3                   
      p(4,3)=p(4,3)+(sisq-cfisq-h)*t4                                   
      p(3,4)=p(3,4)+(cfisq-sisq-h)*t4                                   
      p(4,4)=p(4,4)+(sifi+g)*(t1+t2)+(cfisq+sisq+e)*t3                  
        if(kl.eq.1)then
c         write(6,951)t1,t2,t3,cfisq,sisq,sifi,p(1,1)
951       format('t1,t2,t3,cfisq,sisq,sifi,p11'/1p7e11.3)
        endif
  150 continue                                                          
      do 1 k=1,4                                                        
      kk=(kl-1)*4+k                                                     
      do 1 l=1,4                                                        
      ttpp1=const*p(k,1)*dmu(i)*dp/fi1 
      ttpp2=const*p(k,2)*dmu(i)*dp/fi1 
      ttpp=0.5d0*(ttpp1+ttpp2)
      tp(kk,l)=p(k,l)*dmu(i)*dp/fi1+tp(kk,l)                           
       km=kk-(kl-1)*4
c       if(km.eq.1 .and. l.eq.1)then
c     write(6,952)i,kk,l,dmu(i),p(k,l),tp(kk,l),ttpp
952   format('i,kk,l,dmu,p,tp,ttpp',3i3,1x,1p4e11.3)
c       endif
 1    continue
      do 99 k=1,4                                                       
      kk=(kl-1)*4+k                                                     
      do 99 l=1,4                                                       
      tp(kk,l)=tp(kk,l)/ddmu 
c       if(kk.eq.1 .and. l.eq.1)then
c          write(6,953)kk,l,ddmu,fi1,tp(kk,l)
953        format('kk,l,ddmu,fi1,tp(kk,l)',2i3,1x,1p3e11.3)
c       endif
99     continue                                                         
      id=(kl-1)*4
c      for debug purposes, const=alambd**2/(4*pi**2*qs)
      ppbar=const*(tp(id+1,1)+tp(id+1,2)+tp(id+2,1)+tp(id+2,2))*0.5d0
c     write(6,404)thetd,ppbar, ppbar*pi   
c     do ia=1,4
c        ib=(kl-1)*4+ia
c     write(6,405)(tp(ib,j),j=1,4)
c     average value of the phase function
c     enddo
    7 continue                                                          
c     write(6,*)'leaving hump'
      if(izia.eq.1)stop
c***********************************************************************
c.....format statements.................................................
 23   format(80a1)                                                      
 24   format(10x,80a1)                                                  
 25   format(1p2d18.12,1p2d19.12,0pf6.2)                                
 26   format(1x,1p2d18.12,1p2d19.12,0pf6.2)                             
 404  format(1x,'tp-hump..thetd,ppbar',f8.2,1p2e12.4/)
 405  format(1p4e12.4)                     
 403  format(1x,'t func',i5,1p6e11.3)                                   
c***********************************************************************
      return                                                            
      end                                                               
c***********************************************************************