program hamana2stuff
  implicit none
  real(8)::npix,var1(13),var2(3),x,y,tx,ty,z,kappa,gamma(2)
  character(len=512)::stufffile,shearfile,cnpix
  
  call getarg(1,stufffile)
  call getarg(2,shearfile)
  call getarg(3,cnpix)
  read(cnpix,*)npix

  open(10,file=trim(stufffile))
  open(11,file=trim(shearfile))
  
  do
     read(10,*,end=100)var1(1:13)
     x = var1(2)
     y = var1(3)
     if((x>=0.).and.(y>=0.).and.(x<=npix).and.(y<=npix))then
        read(11,*)var2(1:3)
        kappa    = var2(1)
        gamma(1) = var2(2)
        gamma(2) = var2(3)
        write(100,"(I3,x,2(F10.3,x),F8.4,x,F5.3,5x,2(F5.3,x),F7.2,5x,2(F5.3,x),F7.2,2x,F7.5,2x,F4.1)")&
             & int(var1(1)),var1(2:13)
        call add_shear(var1,gamma,kappa)
        write(*  ,"(I3,x,2(F10.3,x),F8.4,x,F5.3,5x,2(F5.3,x),F7.2,5x,2(F5.3,x),F7.2,2x,F7.5,2x,F4.1)")&
             & int(var1(1)),var1(2:13)
     end if

  end do
100 close(10)
  close(11)

end program hamana2stuff
subroutine add_shear(val,gamma,kappa)
  implicit none
  integer::i,nrow
  real(8)::val(13),posangleB,posangleD,axratB,axratD,radiB,radiD,rotmat(2,2)
  real(8)::rotmatt(2,2),Ellp(2,2),Ell(2,2),matgam(2,2),Ellc(2,2),kappa
  real(8)::lambda(2),x,y,r,rotcos,rotsin,d2r,gamma(2),r1,r2,a,b,c,pi
  character::file*100,cgamm1*10,cgamm2*10
  pi = acos(-1.d0)
  d2r = pi/180.d0
    
  posangleB = val(8)*d2r
  posangleD = val(11)*d2r
  axratB    = val(7)
  axratD    = val(10)
  if(val(6)==0.d0)then
     radiB = 1.d0
  else
     radiB = val(6)
  end if
  if(val(9)==0.d0)then
     radiD = 1.d0
  else
     radiD = val(9)
  end if
  rotmat(1,1) = cos(posangleB)
  rotmat(1,2) = sin(posangleB)
  rotmat(2,1) =-sin(posangleB)
  rotmat(2,2) = cos(posangleB)
  rotmatt(1,1)= cos(posangleB)
  rotmatt(1,2)=-sin(posangleB)
  rotmatt(2,1)= sin(posangleB)
  rotmatt(2,2)= cos(posangleB)
  Ellp(1,1)=radiB/axratB
  Ellp(1,2)=0.d0
  Ellp(2,1)=0.d0
  Ellp(2,2)=radiB*axratB
  Ell = matmul(matmul(rotmatt,Ellp),rotmat)
  matgam(1,1)=1.d0-gamma(1)-kappa
  matgam(1,2)=-gamma(2)
  matgam(2,1)=-gamma(2)
  matgam(2,2)=1.d0+gamma(1)-kappa
  Ellc = matmul(matmul(matgam,Ell),matgam)
  a = 1.d0
  b = -(Ellc(1,1)+Ellc(2,2))
  c = Ellc(1,1)*Ellc(2,2)-Ellc(1,2)*Ellc(2,1)
  lambda(1)=(-b+sqrt(b*b-4.d0*a*c))/2.d0
  lambda(2)=(-b-sqrt(b*b-4.d0*a*c))/2.d0
  x=1.d0
  y=(Ellc(1,1)-lambda(1))/Ellc(1,2)
  r1=sqrt(x**2.d0+y**2.d0)
  y=(Ellc(1,1)-lambda(2))/Ellc(1,2)
  r2=sqrt(x**2.d0+y**2.d0)
  rotcos = 1.d0/r1
  rotsin = (Ellc(1,1)-lambda(1))/Ellc(1,2)/r1
  if(abs(-asin(rotsin)-posangleB).gt.abs(-asin(rotsin)-pi-posangleB))then
     val(8) = (-asin(rotsin)-pi)/d2r
  elseif(abs(-asin(rotsin)-posangleB).gt.abs(-asin(rotsin)+pi-posangleB))then
     val(8) = (-asin(rotsin)+pi)/d2r
  else
     val(8) = -asin(rotsin)/d2r
  end if
  if(val(8).gt. 180.d0)val(8)=val(8)-2.d0*pi/d2r
  if(val(8).lt.-180.d0)val(8)=val(8)+2.d0*pi/d2r
  val(7) = sqrt(minval(lambda)/maxval(lambda))
  if(val(6)/=0.d0)val(6) = Ellp(1,1)*val(7)
  
  rotmat(1,1) = cos(posangleD)
  rotmat(1,2) = sin(posangleD)
  rotmat(2,1) =-sin(posangleD)
  rotmat(2,2) = cos(posangleD)
  rotmatt(1,1)= cos(posangleD)
  rotmatt(1,2)=-sin(posangleD)
  rotmatt(2,1)= sin(posangleD)
  rotmatt(2,2)= cos(posangleD)
  Ellp(1,1)=radiD/axratD
  Ellp(1,2)=0.d0
  Ellp(2,1)=0.d0
  Ellp(2,2)=radiD*axratD
  Ell  = matmul(matmul(rotmatt,Ellp),rotmat)
  Ellc = matmul(matmul(matgam,Ell),matgam)
  a = 1.d0
  b = -(Ellc(1,1)+Ellc(2,2))
  c = Ellc(1,1)*Ellc(2,2)-Ellc(1,2)*Ellc(2,1)
  lambda(1)=(-b+sqrt(b*b-4.d0*a*c))/2.d0
  lambda(2)=(-b-sqrt(b*b-4.d0*a*c))/2.d0
  x=1.d0
  y=(Ellc(1,1)-lambda(1))/Ellc(1,2)
  r1=sqrt(x**2.d0+y**2.d0)
  y=(Ellc(1,1)-lambda(2))/Ellc(1,2)
  r2=sqrt(x**2.d0+y**2.d0)
  rotcos = 1.d0/r1
  rotsin = (Ellc(1,1)-lambda(1))/Ellc(1,2)/r1
  if(abs(-asin(rotsin)-posangleD).gt.abs(-asin(rotsin)-pi-posangleD))then
     val(11) = (-asin(rotsin)-pi)/d2r
  elseif(abs(-asin(rotsin)-posangleD).gt.abs(-asin(rotsin)+pi-posangleD))then
     val(11) = (-asin(rotsin)+pi)/d2r
  else
     val(11) = -asin(rotsin)/d2r
  end if
  if(val(11).gt. 180.d0)val(11)=val(11)-2.d0*pi/d2r
  if(val(11).lt.-180.d0)val(11)=val(11)+2.d0*pi/d2r
  val(10) = sqrt(minval(lambda)/maxval(lambda))
  if(val(9)/=0.d0)val(9) = Ellp(1,1)*val(10)

  return
end subroutine add_shear