********************************** subroutine ang3d(z,nout) ********************************** parameter(ipow=4096,icol=20) implicit double precision(a-h,o-z) dimension z(3,3,2),phi(3),zsz(3),zsqr(3),zsqi(3) dimension zmm(3),zm(3),zmt(3) complex*16 ftest,fi,fb3,fb2,fb1 real*4 ans(ipow,icol) common/result/ans zpi=8.*atan(1.d0) c... calculate amplitude and phase of all components, c... real and imaginary parts of squares of all components do 10 i=1,3 zsz(i)=z(i,3,1)*z(i,3,1)+z(i,3,2)*z(i,3,2) phi(i)=datan2(z(i,3,2),z(i,3,1)) zsqr(i)=z(i,3,1)**2-z(i,3,2)**2 10 zsqi(i)=2.d0*z(i,3,1)*z(i,3,2) c... calculate real and imaginary part of (z1**2+z2**2+z3**2) c... and amplitude and phase zzr=zsqr(1)+zsqr(2)+zsqr(3) zzi=zsqi(1)+zsqi(2)+zsqi(3) * zzsq=zzr*zzr+zzi*zzi * zzsq=sqrt(zzsq) thd=.5*atan2(zzi,zzr) r1=zsz(1)*cos(thd-phi(1))**2+zsz(2)*cos(thd-phi(2))**2+ & zsz(3)*cos(thd-phi(3))**2 if((zsz(1)+zsz(2)+zsz(3)).gt.3.d0*r1) thd=thd+zpi/4. c... calculate polarization basis vectors in real 3-D space cs=cos(thd) sn=sin(thd) zmm(1)=z(1,3,1)*cs+z(1,3,2)*sn zmm(2)=z(2,3,1)*cs+z(2,3,2)*sn zmm(3)=z(3,3,1)*cs+z(3,3,2)*sn bzmm=sqrt(zmm(1)*zmm(1)+zmm(2)*zmm(2)+zmm(3)*zmm(3)) do 20 i=1,3 20 zmm(i)=zmm(i)/bzmm cs=cos(thd-zpi/4.) sn=sin(thd-zpi/4.) zm(1)=z(1,3,1)*cs+z(1,3,2)*sn zm(2)=z(2,3,1)*cs+z(2,3,2)*sn zm(3)=z(3,3,1)*cs+z(3,3,2)*sn bzm=sqrt(zm(1)*zm(1)+zm(2)*zm(2)+zm(3)*zm(3)) do 30 i=1,3 30 zm(i)=zm(i)/bzm c... major axis should have z < 0 if(zmm(1).gt.-1.e-09)then zmm(1)=-zmm(1) zmm(2)=-zmm(2) zmm(3)=-zmm(3) endif c... minor axis should have z > 0 if(zm(1).lt.1.e-09)then zm(1)=-zm(1) zm(2)=-zm(2) zm(3)=-zm(3) endif c... right hand rule for normal vector on ellipse zmt(1)=zmm(2)*zm(3)-zmm(3)*zm(2) zmt(2)=zmm(3)*zm(1)-zmm(1)*zm(3) zmt(3)=zmm(1)*zm(2)-zmm(2)*zm(1) c.................................... c... test for re(arctan(z3/z2)) fi=cmplx(0.,1.) ftest=cmplx(z(3,3,1),z(3,3,2))/cmplx(z(2,3,1),z(2,3,2)) ftest=(cdlog(1.+fi*ftest)-cdlog(1-fi*ftest))/2./fi frtest=real(ftest)*360./zpi fb1=cmplx(z(1,3,1),z(1,3,2)) fb2=cmplx(z(2,3,1),z(2,3,2)) fb3=cmplx(z(3,3,1),z(3,3,2)) ftest=fb3/fb2 ftest=atan(cdabs(fb3)/cdabs(fb2)) c ftest=(cdlog(1.+fi*ftest)-cdlog(1-fi*ftest))/2./fi ... arctan c ftest=-fi*cdlog(ftest+fi*cdsqrt(1-ftest*ftest)) ... arccos fitest=real(ftest) c ftest=-fi*cdlog(fi*ftest+cdsqrt(1-ftest*ftest)) ... arcsin fltest=dimag(ftest) fltest=atan2(fitest,fltest)*360/zpi c print 1800,' freq, re(arctan(z3/z2)), ... ',ans(nout,1),frtest, c & fitest,fltest c1800 format(a30,e15.5,3f10.6) c... end of test c.................................... c... calculation of tilt, use minor axis for calculation if c... major axis is along vert. axis c... sign is pos. if looking downward the major axis and ellipse c... is tilted to the right if(abs(zmm(2)).lt.1.e-03.and.abs(zmm(3)).lt.1.e-03)then ans(nout,16)=360.d0/zpi*atan2(zm(3),zm(2)) else fnom=(zmm(3)*zmt(2)-zmm(2)*zmt(3))/ & sqrt(zmm(2)*zmm(2)+zmm(3)*zmm(3)) ans(nout,16)=dsign(1.d0,(zm(1)*zmt(1)))*360.d0/zpi*acos(fnom) endif c... calculation of the two principal angles of major axis in 3-d space ans(nout,17)=360.d0/zpi*atan2(zmm(3),zmm(2)) if(abs(zmm(3)).lt.1.e-04)ans(nout,17)=0. fnom=sqrt(zmm(3)*zmm(3)+zmm(2)*zmm(2)) ans(nout,18)=360.d0/zpi*atan2(fnom,-zmm(1)) c... calculation of ellipticity and noise-to-signal ratio (paulssen c... et al.) ans(nout,19)=1.d0-bzm/bzmm ans(nout,20)=(1.d0-ans(nout,7))/(3.d0*ans(nout,7)-1.d0) return end ***************************************************** subroutine eptis(iwave,nout) ***************************************************** c... evaluation of ellipticities c... unwrapping of angles c... this code changes the content coming out of POLAR! c... the new polarization parameters will be as following: c output c all results in common/result/ans(1.....nout,1-20) c ans(..,1)=frequency (in Hz) c ans(..,2)=amplitude of component 1, z c ans(..,3)=amplitude of component 2, r c ans(..,4)=amplitude of component 3, t c ans(..,5)=smallest sing val c ans(..,6)=middle sing val c ans(..,7)=largest sing val c ans(..,8)=horiz azimuth c ans(..,9)= and its error c ans(..,10)=sign of horizontal phase (+prograde,-retrograde) c ans(..,11)= error of eps-HH (eval. from error of phi-HH) c ans(..,12)=vertical azimuth (same as ans(..,18)) c ans(..,13)= and its error c ans(..,14)= sign of vertical phase (+prograde,-retrograde) c ans(..,15)= error of eps-VH (eval. from error of phi-VH) c...ans(..,16)= tilt, theta-3T c...ans(..,17)= hor. proj. of major axis, theta-3H c...ans(..,18)= vert. angle of major axis, theta-3V c...ans(..,19)= ellipticity (1-b/a) c...ans(..,20)= used by polart for wrapping ans(17) parameter(ipow=4096,icol=20) double precision rat,phihh,phivh,zpi double precision sqhh,sqvh common/result/ans(ipow,icol) common/resule/epshh(ipow),epsvh(ipow) yminh0=-180. alim=90. zpi=4.*atan(1.d0) do 50 i=1,nout rat=ans(i,3)/ans(i,4) phihh=zpi/360.*ans(i,10) phivh=zpi/360.*ans(i,14) ans(i,11)=zpi/360.*ans(i,11) ans(i,15)=zpi/360.*ans(i,15) sqhh=sqrt((rat*rat+1)**2-4.*sin(phihh)*sin(phihh)*rat*rat) epshh(i)=2*sqhh/(rat*rat+1+sqhh) ans(i,11)=ans(i,11)*4.*rat*rat*(rat*rat+1)*sin(2*phihh)/ & (rat*rat+1+sqhh)**2/sqhh rat=ans(i,2)/sqrt(ans(i,3)*ans(i,3)+ans(i,4)*ans(i,4)) sqvh=sqrt((rat*rat+1)**2-4.*sin(phivh)*sin(phivh)*rat*rat) epsvh(i)=2*sqvh/(rat*rat+1+sqvh) ans(i,15)=ans(i,15)*4.*rat*rat*(rat*rat+1)*sin(2*phivh)/ & (rat*rat+1+sqvh)**2/sqvh ans(i,10)=sign(1.,ans(i,10)) ans(i,14)=sign(1.,ans(i,14)) ans(i,20)=0. c... wrapping of angles if(iwave.gt.0)then goto (51,52,53) iwave 51 continue if(abs(ans(i,8)).gt.abs(yminh0/2)) & ans(i,8)=ans(i,8)+sign(1,ans(i,8))*yminh0 if(abs(ans(i,17)).gt.abs(yminh0/2))then ans(i,17)=ans(i,17)+sign(1,ans(i,17))*yminh0 ans(i,20)=1. endif goto 50 52 continue if(ans(i,17).lt.0.)then ans(i,17)=ans(i,17)+180. ans(i,20)=1. endif if(ans(i,8).lt.0.) & ans(i,8)=ans(i,8)+180. goto 50 53 continue if(i.lt.2)then if(abs(ans(i,17)).gt.180.) & ans(i,17)=ans(i,17)-180.*sign(1.,ans(i,17)) if(abs(ans(i,8)).gt.180.) & ans(i,8)=ans(i,8)-180.*sign(1.,ans(i,8)) else if(abs(ans(i,17)-ans(i-1,17)).gt.alim)then ans(i,17)=ans(i,17)-180.* & sign(1.,ans(i,17)-ans(i-1,17)) ans(i,20)=1. endif if(abs(ans(i,8)-ans(i-1,8)).gt.alim) & ans(i,8)=ans(i,8)-180.* & sign(1.,ans(i,8)-ans(i-1,8)) endif goto 50 endif 50 continue if(iwave.eq.3)then do 60 i=2,nout if(abs(ans(i,17)-ans(i-1,17)).gt.alim)then ans(i,17)=ans(i,17)-180.* & sign(1.,ans(i,17)-ans(i-1,17)) if(ans(i,20).gt.0.9)then ans(i,20)=0. else ans(i,20)=1. endif endif if(abs(ans(i,8)-ans(i-1,8)).gt.alim) & ans(i,8)=ans(i,8)-180.* & sign(1.,ans(i,8)-ans(i-1,8)) 60 continue endif return end