!     Last change:  ER    9 Nov 2006    1:55 pm
      subroutine polfa (thet,thet1,thet2,e0x,e0y,e0z,pk_01,pk_02,pk_12,pk_21,cos_01_21)
      use intensity_data, only: pi,iunit5,iunit6

      del=(cos(2.*thet2)-cos(2.*thet1)*cos(2.*thet))/(sin(2.*thet1)*sin(2.*thet))
      if(abs(del)>1.)del=sign(1.,del)
      del=acos(del)
      
      hmue=asin(sin(2.*thet1)*sin(del))
      if(del.gt.pi/2.)hmue=pi-hmue

if(abs(cos(2.*thet1)/cos(hmue))>1.)then
      thx=acos(sign(1.,cos(2.*thet1)/cos(hmue)))
else
      thx=acos(cos(2.*thet1)/cos(hmue))
endif
 !!!  s1
      xkco=cos(2.*thet)
      ykco=sin(2.*thet)
      zkco=0.
 !!!  s2
      xkop=cos(2.*thet1)
      ykop=sin(thx)*cos(hmue)
      zkop=sin(hmue)

 !!!  s2 x E0
      zwix=ykop*e0z-zkop*e0y
      zwiy=zkop*e0x-xkop*e0z
      zwiz=xkop*e0y-ykop*e0x
 !!!  - s2  x (s2 x E0)
      e1x=-(ykop*zwiz-zkop*zwiy)
      e1y=-(zkop*zwix-xkop*zwiz)
      e1z=-(xkop*zwiy-ykop*zwix)
      e1l=e1x**2+e1y**2+e1z**2
      pk_02=SQRT(e1l)
 !!! unit vector E2
      e1x=e1x/sqrt(e1l)
      e1y=e1y/sqrt(e1l)
      e1z=e1z/sqrt(e1l)
 !!!  s1 x E2
      zwix=ykco*e1z-zkco*e1y
      zwiy=zkco*e1x-xkco*e1z
      zwiz=xkco*e1y-ykco*e1x
 !!!  - s1  x (s1 x E2)
      e2x=-(ykco*zwiz-zkco*zwiy)
      e2y=-(zkco*zwix-xkco*zwiz)
      e2z=-(xkco*zwiy-ykco*zwix)
      e2l=e2x**2+e2y**2+e2z**2
      pk_21=SQRT(e2l)

 !!!  s1 x E0
      zwix=ykco*e0z-zkco*e0y
      zwiy=zkco*e0x-xkco*e0z
      zwiz=xkco*e0y-ykco*e0x
 !!!  - s1  x (s1 x E0)
      e1x=-(ykco*zwiz-zkco*zwiy)
      e1y=-(zkco*zwix-xkco*zwiz)
      e1z=-(xkco*zwiy-ykco*zwix)
      e1l=e1x**2+e1y**2+e1z**2
      pk_01=SQRT(e1l)
 !!!  cos des Winkels zwischen E1 und E-Umweg 
      cos_01_21=(e1x*e2x+e1y*e2y+e1z*e2z)/sqrt(e1l*e2l)

  !!! unit vector E1
      e1x=e1x/sqrt(e1l)
      e1y=e1y/sqrt(e1l)
      e1z=e1z/sqrt(e1l)
 !!!  s2 x E1
      zwix=ykop*e1z-zkop*e1y
      zwiy=zkop*e1x-xkop*e1z
      zwiz=xkop*e1y-ykop*e1x
 !!!  - s2  x (s2 x E1)
      e2x=-(ykop*zwiz-zkop*zwiy)
      e2y=-(zkop*zwix-xkop*zwiz)
      e2z=-(xkop*zwiy-ykop*zwix)
      e2l=e2x**2+e2y**2+e2z**2
      pk_12=SQRT(e2l)

      return
      end
      subroutine upsi(u,v,w,ub)
      use intensity_data, only: h,k,l,psid

implicit none
real    :: c,p,ps,r1,r2,r3,r4,r5,r6,r7,r8,r9
real    :: t1,t2,t3,t4,t5,t6,tl,u,v,w,uvw1,uvw2,uvw3,uvwl 
real    :: ub(9)

      t1=ub(1)*h+ub(2)*k+ub(3)*l
      t2=ub(4)*h+ub(5)*k+ub(6)*l
      t3=ub(7)*h+ub(8)*k+ub(9)*l
      tl=sqrt(t1*t1+t2*t2+t3*t3)
      t1=t1/tl
      t2=t2/tl
      t3=t3/tl
      if(t1.eq.0.and.t2.eq.0.)then
      p=0.
      else
      p=mod(6.28318-atan2(t1,t2),6.28318)
      endif
      c=atan2(t3,sqrt(t1*t1+t2*t2))
      t4=ub(1)*u+ub(2)*v+ub(3)*w
      t5=ub(4)*u+ub(5)*v+ub(6)*w
      t6=ub(7)*u+ub(8)*v+ub(9)*w
      tl=sqrt(t4*t4+t5*t5+t6*t6)
      t4=t4/tl
      t5=t5/tl
      t6=t6/tl
      r1=cos(p)
      r2=sin(p)
      r3=0.
      r4=-cos(c)*sin(p)
      r5=cos(c)*cos(p)
      r6=sin(c)
      r7=sin(c)*sin(p)
      r8=-sin(c)*cos(p)
      r9=cos(c)
      uvw1=r1*t4+r2*t5+r3*t6
      uvw2=r4*t4+r5*t5+r6*t6
      uvw3=r7*t4+r8*t5+r9*t6
      uvwl=sqrt(uvw1*uvw1+uvw2*uvw2+uvw3*uvw3)
      ps=acos(uvw1/uvwl)
      psid=ps*57.2958*sign(1.,uvw3)
      return
      end