!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function rotsou2_geodip(Ccig) result(Scig)

 complex				:: Ccig(:,:,:)
 complex, allocatable			:: Scig(:,:,:)
 integer				:: ix, iz, ih, stat, ikhrot,k
 real					:: tgama,talfa,tteta_s,gwrot,fwrot,khrot

 allocate(Scig(kz%n,kh%n,kx%n));Scig=0.

k=0
 do ix=1,kx%n
    do iz=1,kz%n
       if(kz%all(iz).eq.0.) then
          if(kx%all(ix).lt.0.) then
             talfa=1.E+37
          else if(kx%all(ix).gt.0.) then
             talfa=-1.E+37
          else
             talfa=-kx%d/kz%d
          end if
       else
          talfa=-kx%all(ix)/kz%all(iz)
       end if

       do ih=1,kh%n
          if(kz%all(iz).eq.0.) then
             if(kx%all(ix).lt.0.) then
                khrot=kh%all(ih)
             else if(kx%all(ix).gt.0.) then
                khrot=-kh%all(ih)
             end if
          else
             tteta_s=-kh%all(ih)/kz%all(iz)
          end if

	  if (abs(1.+talfa*tteta_s).lt.(1.E-37)) then
             if(kx%all(ix).lt.0.) then
                khrot=kh%all(ih)
             else if(kx%all(ix).gt.0.) then
                khrot=-kh%all(ih)
             end if
	  else
             tgama=(tteta_s-talfa)/(1.+talfa*tteta_s)
	     khrot=-kz%all(iz)*tgama
	  end if

	  ikhrot=int((khrot-kh%o)/kh%d)+1

	  fwrot=(khrot-kh%o)/kh%d-ikhrot-1
	  gwrot=1.-fwrot

  	  if (ikhrot.ge.1 .and. ikhrot.lt.kh%n) then
	     k=k+1
  	     Scig(iz,ih,ix)=gwrot*Ccig(iz,ikhrot,ix)+fwrot*Ccig(iz,ikhrot+1,ix)+Scig(iz,ih,ix)
	  end if

       end do
    end do
 end do
 write(0,*) k

 end function rotsou2_geodip

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 function rotsou1_geodip(Ccig) result(Scig)

 complex				:: Ccig(:,:,:)
 complex, allocatable			:: Scig(:,:,:)
 integer				:: ix, iz, ih, stat, ikhrot1,ikhrot2,k
 real					:: tgama,talfa,tteta_s,gwrot,fwrot,khrot
 real                                   :: lkh, lkx, lkz

 allocate(Scig(kz%n,kh%n,kx%n));Scig=0.

k=0

 do ix=1,kx%n
    lkx = kx%all(ix)
    if (lkx > pi/dh) lkx = lkx-2*pi/dh
 
    do iz=1,kz%n
       lkz = kz%all(iz)
       if (lkz > pi/dh) lkz = lkz-2*pi/dh

       if(lkz.eq.0.) then
          if     (lkx.lt.0.) then; talfa=1.E+37
          else if(lkx.gt.0.) then; talfa=-1.E+37
          else                   ; talfa=-kx%d/kz%d
          end if
       else
          talfa=-lkx/lkz
       end if

       do ih=1,kh%n
          lkh = kh%all(ih)
          if (lkh > pi/dh) lkh = lkh-2*pi/dh

          if(lkz.eq.0.) then
             if     (lkh.lt.0.) then; tteta_s=1.E+37
             else if(lkh.gt.0.) then; tteta_s=-1.E+37
             else		    ; tteta_s=kh%d/kz%d
             end if
          else
             tteta_s=-lkh/lkz
          end if

	  if (abs(1.+talfa*tteta_s).lt.(1.E-37)) then
	     if(tteta_s-talfa.lt.0.) then
	       tgama=-1.E+37
	     else if(tteta_s-talfa.gt.0.) then
	       tgama=1.E+37
	     else
	       tgama=1.
	     end if
	  else
             tgama=(tteta_s-talfa)/(1.+talfa*tteta_s)
	  end if

	  khrot=-lkz*tgama
	  ikhrot1=int(khrot/kh%d-pi/dh)+1
!if (khrot < 0) write(0,*) khrot

  	  if (ikhrot1.ge.1 .and. ikhrot1.lt.kh%n) then
	     ikhrot2=mod(ikhrot1,kh%n)+1
	     fwrot=khrot/kh%d-ikhrot1-1
	     gwrot=1.-fwrot

	     k=k+1
  	     Scig(iz,ih,ix)=gwrot*Ccig(iz,ikhrot1,ix)+fwrot*Ccig(iz,ikhrot2,ix)+Scig(iz,ih,ix)
	  end if

       end do
    end do
 end do
 write(0,*) k

 end function rotsou1_geodip