SUBROUTINE PRINAXIS(MXNAT,NAT,ICOMP,IXYZ,A1,A2) C Arguments: INTEGER MXNAT,NAT INTEGER*2 IXYZ(MXNAT,3),ICOMP(MXNAT,3) REAL A1(3),A2(3) C C Local variables: INTEGER LWORK PARAMETER(LWORK=12) INTEGER & I,INFO, & J,J1,J2,J3,JA,JD, & K,LDA,LDVL,LDVR,N REAL DENOM,PI,ZNORM,ZNORM2 REAL & A(3,3),A3(3),VL(3,3),VR(3,3),WI(3),WORK(LWORK),WR(3), & ZI(3,3),ZR(3,3) DOUBLE PRECISION DSUM DOUBLE PRECISION CM(3) CHARACTER CMSGNM*70 C*********************************************************************** C Subroutine to determine principal axes of moment of inertia tensor C of target. C Target is assumed to be homogeneous. C Given: C C MXNAT =limit on largest allowed value of NAT=number of dipoles C NAT =number of dipoles in target C IXYZ(J,1-3) =x,y,z location of dipole J on lattice C ICOMP(J,1-3)=composition for dipole J; x,y,z directions C C Returns: C C A1(1-3) =principal axis A1 in target frame (normalized) C A2(1-3) =principal axis A2 in target frame (normalized) C C where A1 = principal axis with largest moment of inertia C A2 = principal axis with second-largest moment of inertia C C B.T. Draine, Princeton Univ. Observatory, 96.01.26 C History: C 96.01.26 (BTD) extracted from routine TARBLOCKS C IX,IY,IZ replaced by IXYZ C 96.02.23 (BTD) corrected typo in 2 format statements C add code computing axis a3 and printing out C eigenvalue. C 96.11.01 (PJF) modified to use LAPACK code rather than EIGV C and associated routines from NSWC package. C 96.11.21 (BTD) modified to remove WRITE(0 and replace with C calls to WRIMSG C 98.12.15 (BTD) Appears to be a problem with LAPACK routine SGEEV C when compiled using g77 under Linux. C Add trap to catch this problem if it occurs, C print out warning, and terminate execution. C End history. C C Copyright (C) 1993,1994,1995,1996,1998 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) C C Compute center of mass (assuming uniform density). C DO JD=1,3 CM(JD)=0. ENDDO DO JA=1,NAT DO JD=1,3 CM(JD)=CM(JD)+IXYZ(JA,JD) ENDDO ENDDO DO JD=1,3 CM(JD)=CM(JD)/NAT ENDDO C C Compute moment of inertia tensor C DO J=1,3 DO K=1,3 A(J,K)=0. ENDDO ENDDO DSUM=0. DO K=1,3 DO JA=1,NAT DSUM=DSUM+(IXYZ(JA,K)-CM(K))**2 ENDDO ENDDO DO J=1,3 A(J,J)=DSUM+NAT/6. ENDDO DO J=1,3 DO K=1,3 DSUM=0. DO JA=1,NAT DSUM=DSUM+(IXYZ(JA,J)-CM(J))*(IXYZ(JA,K)-CM(K)) ENDDO A(J,K)=A(J,K)-DSUM ENDDO ENDDO C C Normalize moment of inertia tensor to units of 0.4*mass*a_eff**2 C DENOM=0.4*NAT*(.75*NAT/PI)**(2./3.) DO 3500 K=1,3 DO 3400 J=1,3 A(J,K)=A(J,K)/DENOM 3400 CONTINUE 3500 CONTINUE c*** diagnostic c write(0,*)'in prinaxis:' c write(0,*)'center of mass CM=',CM c write(0,*)'--- normalized moment of inertia tensor A ---' c write(0,*)(a(1,k),k=1,3) c write(0,*)(a(2,k),k=1,3) c write(0,*)(a(3,k),k=1,3) c*** C C Compute eigenvalues and eigenvectors of A(J,K) C N=3 LDA=3 LDVR=3 LDVL=3 C calling sequence of SGEEV from "lapack.f" C CALL SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, C $ LDVR, WORK, LWORK, INFO ) C C in fact only right eignevectors are needed C CALL SGEEV( 'V', 'V', N, a, LDA, WR, WI, VL, LDVL, VR, $ LDVR, WORK, LWORK, INFO ) C 98.12.06 BTD It appears that SGEEV or subsidiary LAPACK routines C are not executing properly under Linux g77. This C problem has not been noticed previously, and certainly C did not appear with the Solaris f77 compiler. Perhaps C one of the LAPACK routines is failing to SAVE some C variable between calls. In any event, we have C introduced this trap to terminate execution in case C of failure of SGEEV. C IF(INFO.NE.0)THEN WRITE(CMSGNM,FMT='(A,I4)') & 'Error in PRINAXIS: returned from SGEEV with INFO=',INFO WRITE(CMSGNM,FMT='(A)') & 'SGEEV is part of lapack.f' WRITE(CMSGNM,FMT='(A)') & 'LAPACK routines may not work under Linux g77' WRITE(CMSGNM,FMT='(A)') & 'workaround by specifying target axis by hand' CALL ERRMSG('FATAL','PRINAXIS',' failure in SGEEV') ENDIF C c VR (output) REAL array, dimension (LDVR,N) c If JOBVR = 'V', the right eigenvectors v(j) are stored one c after another in the columns of VR, in the same order c as their eigenvalues. c If JOBVR = 'N', VR is not referenced. c If the j-th eigenvalue is real, then v(j) = VR(:,j), c the j-th column of VR. c If the j-th and (j+1)-st eigenvalues form a complex c conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and c v(j+1) = VR(:,j) = i*VR(:,j+1). C Copy to zr and zi arrays. This is for historical reasons C because of previous use of "NSWC Library of mathematical subroutines". C (code eigv). There is a slight differnce between "eigv" C and "sgeev" drivers in that eigenvectors in "eigv" are returned C as "complex" and "real" part but "sgeev" return just one matrix. C Now, for real eigenvalues the complex part of eignevector should C be 0. Thus, there should not be a difference between "eigv" C and "sgeev" driver. C do i=1,3 do j=1,3 zr(i,j)=vr(i,j) zi(i,j)=0. enddo enddo C C Normalize the eigenvectors: C DO 4000 J=1,3 ZNORM2=0. DO 3800 K=1,3 ZNORM2=ZNORM2+ZR(K,J)**2+ZI(K,J)**2 3800 CONTINUE C C Make first component of each eigenvector positive C ZNORM=SQRT(ZNORM2) IF(ZR(1,J).LT.0.)THEN ZNORM=-ZNORM ELSEIF(ZR(1,J).EQ.0.)THEN IF(ZR(2,J).LT.0.)ZNORM=-ZNORM ENDIF DO 3900 K=1,3 ZR(K,J)=ZR(K,J)/ZNORM ZI(K,J)=ZI(K,J)/ZNORM 3900 CONTINUE 4000 CONTINUE C C Order the eigenvectors by eigenvalue, largest first: C J1=1 IF(WR(2).GT.WR(1))J1=2 IF(WR(3).GT.WR(J1))J1=3 IF(J1.EQ.1)THEN J2=2 IF(WR(3).GT.WR(2))J2=3 ELSEIF(J1.EQ.2)THEN J2=1 IF(WR(3).GT.WR(1))J2=3 ELSEIF(J1.EQ.3)THEN J2=1 IF(WR(2).GT.WR(1))J2=2 ENDIF DO 4100 K=1,3 A1(K)=ZR(K,J1) A2(K)=ZR(K,J2) 4100 CONTINUE C C Set J3=3 if J1+J2=1+2=3 C 2 =1+3=4 C 1 =2+3=5 C J3=6-J1-J2 C C Having fixed eigenvectors a_1 and a_2, compute a_3=a_1xa_2 C A3(1)=A1(2)*A2(3)-A1(3)*A2(2) A3(2)=A1(3)*A2(1)-A1(1)*A2(3) A3(3)=A1(1)*A2(2)-A1(2)*A2(1) C C*** enable following lines to print out eigenvalues WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_1/.4Ma_eff^2 =',WR(J1) CALL WRIMSG('PRINAXIS',CMSGNM) WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_2/.4Ma_eff^2 =',WR(J2) CALL WRIMSG('PRINAXIS',CMSGNM) WRITE(CMSGNM,FMT='(A,0PF8.5)')'I_3/.4Ma_eff^2 =',WR(J3) CALL WRIMSG('PRINAXIS',CMSGNM) WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')'axis a_1 = (',A1,')' CALL WRIMSG('PRINAXIS',CMSGNM) WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')' a_2 = (',A2,')' CALL WRIMSG('PRINAXIS',CMSGNM) WRITE(CMSGNM,FMT='(A,0P3F9.5,A)')' a_3 = (',A3,')' CALL WRIMSG('PRINAXIS',CMSGNM) C*** RETURN END SUBROUTINE ROT2(A,THETA,RM) C*********************************************************************** C This subroutine finds the (3 by 3) rotation matrix for the C given axis (vector) --- a and the rotation angle --- theta. C On input: C a(3) --- axis of rotation (vector, any length) C theta --- angle of rotation (in radians) around vector a C On output: C rm(3,3) --- rotation matrix such that new vector v2 is obtained C from the old vector v1 by the transformation: C v2 = rm \dot v1 C (If more rotations are needed: generate matrices rm with the C help of this subroutine and (matrix) multiply them before C performing v2 = rm \dot v1.) C C Reference: C Leubner, C., 1977, Coordinate-free rotation operator, C Am. J. Phys. 47, 727---729. C C History: C Written by PJF, May 1989, for use in the rotation-averaged C DDA calculations. C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C pre-calculate some stuff: C .. Scalar Arguments .. REAL THETA C .. C .. Array Arguments .. REAL A(3), RM(3,3) C .. C .. Local Scalars .. REAL ANORM, CT, ST INTEGER I, J C .. C .. Local Arrays .. REAL AHAT(3) C .. C .. External Subroutines .. EXTERNAL XNORM3 C .. C .. Intrinsic Functions .. INTRINSIC COS, SIN C .. CT = COS(THETA) ST = SIN(THETA) CALL XNORM3(A,ANORM) AHAT(1) = A(1)/ANORM AHAT(2) = A(2)/ANORM AHAT(3) = A(3)/ANORM C cos(\theta) {\bf 1} term: RM(1,1) = CT RM(2,2) = CT RM(3,3) = CT C Skew-term a \times {\bf 1}: RM(1,2) = -ST*AHAT(3) RM(1,3) = ST*AHAT(2) RM(2,1) = ST*AHAT(3) RM(2,3) = -ST*AHAT(1) RM(3,1) = -ST*AHAT(2) RM(3,2) = ST*AHAT(1) C aa-dyadic (outer-product): DO 20 I = 1, 3 DO 10 J = 1, 3 RM(I,J) = RM(I,J) + (1.-CT)*AHAT(J)*AHAT(I) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE XNORM3(A,XN) C auxiliary for "rotation" routines C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C .. Scalar Arguments .. REAL XN C .. C .. Array Arguments .. REAL A(3) C .. C .. Intrinsic Functions .. INTRINSIC SQRT C .. XN = SQRT(A(1)**2+A(2)**2+A(3)**2) RETURN END SUBROUTINE MULT3(A,B,C) C*********************************************************************** C Purpose: C Given: A(I,J)= 3x3 matrix C B(J,K)= 3x3 matrix C Returns:C(I,K)= matrix product of A and B C This is auxilary for "rotation" routines C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. *********************************************************************** C Arguments: REAL A(3,3),B(3,3),C(3,3) C Local variables: REAL S INTEGER I,J,K C DO 30 J = 1, 3 DO 20 I = 1, 3 S = 0. DO 10 K = 1, 3 S = S + A(I,K)*B(K,J) 10 CONTINUE C(I,J) = S 20 CONTINUE 30 CONTINUE RETURN END SUBROUTINE PROD3(RM,VIN,VOUT) C Given: RM(I,J)=3x3 matrix C VIN(J)=3-vector C Returns:VOUT(I)=RM*VIN C C Auxiliary for "rotation" routines C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C Arguments; REAL RM(3,3), VIN(3), VOUT(3) C Local variables: INTEGER I,J C DO 20 I = 1, 3 VOUT(I) = 0. DO 10 J = 1, 3 VOUT(I) = VOUT(I) + RM(I,J)*VIN(J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE PROD3C(RM,CVIN,CVOUT) C Arguments: REAL RM(3,3) COMPLEX CVIN(3),CVOUT(3) C Local variables: INTEGER I,J C Given: C RM=3x3 rotation matrix C CVIN=complex 3-vector C Computes: C CVOUT=RM \dot CVIN C C Copyright (C) 1993, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** DO 20 I=1,3 CVOUT(I)=(0.,0.) DO 10 J=1,3 CVOUT(I)=CVOUT(I)+RM(I,J)*CVIN(J) 10 CONTINUE 20 CONTINUE RETURN END SUBROUTINE PROD3V(RM,VIN,VOUT,NV,NVMAX) C Arguments: INTEGER NV,NVMAX REAL RM(3,3),VIN(3,NVMAX),VOUT(3,NVMAX) C Local variables: INTEGER I,J,N C Given: C RM = 3x3 rotation matrix C VIN = NV 3-vectors C Computes: C VOUT=RM \dot VIN for each of the NV 3-vectors C C Copyright (C) 1993,1996 B.T. Draine and P.J. Flatau C History C 96.08.09 (BTD) trivial modification to remove obsolescent construct. C end history C This code is covered by the GNU General Public License. C*********************************************************************** DO 20 I=1,3 DO 10 N=1,NV VOUT(I,N)=0. 10 CONTINUE 20 CONTINUE C DO 50 I=1,3 DO 40 J=1,3 DO 30 N=1,NV VOUT(I,N)=VOUT(I,N)+RM(I,J)*VIN(J,N) 30 CONTINUE 40 CONTINUE 50 CONTINUE RETURN END SUBROUTINE ROTATE(A1,A2,AK1,CXE01,CXE02,BETA,THETA,PHI, & EN0R,CXE01R,CXE02R,MXSCA,NSCAT,ENSC,EM1,EM2, & AKSR,EM1R,EM2R) C Arguments: INTEGER MXSCA,NSCAT REAL AK1, BETA, THETA, PHI REAL A1(3), A2(3), AKSR(3,MXSCA), & EM1(3,MXSCA), EM1R(3,MXSCA), EM2(3,MXSCA), & EM2R(3,MXSCA), EN0R(3), ENSC(3,MXSCA) COMPLEX CXE01(3),CXE02(3),CXE01R(3),CXE02R(3) C Local variables INTEGER I REAL BETA0, PHI0, SINTHE, THETA0 REAL RM1(3,3),RM2(3,3),RM3(3,3),VEC(3) C*********************************************************************** C Given: C A1(1-3)=axis 1 of (original, unrotated) target C A2(1-3)=axis 2 of (original, unrotated) target C AK1=magnitude of k vector C CXE01(1-3)=original incident polarization vector 1 C CXE02(1-3)=original incident polarization vector 2 C BETA=angle (radians) through which target is to be rotated C around target axis A1 C THETA=angle (radians)which target axis A1 is to make with C Lab x-axis C PHI=angle (radians) between Lab x,A1 plane and Lab x,y plane C ENSC(1-3,1-NSCAT)=scattering directions in Lab Frame C EM1(1-3,1-NSCAT)=scattering polarization vector 1 in Lab Frame C EM2(1-3,1-NSCAT)=scattering polarization vector 2 in Lab Frame C C Returns: C EN0R(1-3)=incident propagation vector in Target Frame C CXE01R(1-3)=incident polariz. vector 1 in Target Frame C CXE02R(1-3)=incident polariz. vector 2 in Target Frame C AKSR(1-3,1-NSCAT)=scattering k vectors in Target Frame C EM1R(1-3,1-NSCAT)=scattering pol. vector 1 in Target Frame C EM2R(1-3,1-NSCAT)=scattering pol. vector 2 in Target Frame C C Purpose: C In the Lab Frame, we hold the propagation direction (x-axis) C and incident polarization vectors (CXE01 and CXE02) fixed, C and rotate target to an orientation specified by the three C angles BETA,THETA,PHI. C However, computations in the main program DDSCAT are carried out C in the Target Frame, in which the target is held fixed and the C propagation vectors and polarization vectors are rotated. C Hence, this routine computes the incident propagation vector, C scattering propagation vectors, and associated polarization C vectors in the Target Frame. C C B. T. Draine, Princeton Univ. Obs., 89.11.20 C History: C 91.04.22 (BTD) Eliminated vector A1R and removed superfluous line C computing A1R C 96.11.03 (BTD) Corrected comments C end history C C Copyright (C) 1996, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C** Determine the orientation angles for original unrotated target: THETA0=ACOS(A1(1)) SINTHE=SIN(THETA0) IF(SINTHE.NE.0.)THEN PHI0=ACOS(A1(2)/SIN(THETA0)) ELSE C** If sin(THETA0)=0, we will arbitrarily assume PHI0=0 PHI0=0. ENDIF IF(SINTHE.NE.0.)THEN BETA0=ACOS(-A2(1)/SINTHE) ELSE BETA0=ACOS(A2(2)) ENDIF C C**** First rotate the target through angle PHI-PHI0 around x axis C (i.e., rotate the lab through angle PHI0-PHI around x axis) C** obtain rotation matrix: VEC(1)=1. VEC(2)=0. VEC(3)=0. CALL ROT2(VEC,PHI0-PHI,RM1) C C**** Now rotate the target through angle THETA-THETA0 around axis C VEC = xaxis cross a1 C (i.e., rotate the lab through angle THETA0-THETA around C xaxis cross a1 ) C** first compute rotation axis VEC(1)=0. VEC(2)=-A1(3) VEC(3)=A1(2) C*** For special case where a1=xaxis, we want rotation axis=a1 cross a2 IF((VEC(2)**2+VEC(3)**2).LT.1.E-10)THEN VEC(2)=A1(3)*A2(1)-A1(1)*A2(3) VEC(3)=A1(1)*A2(2)-A1(2)*A2(1) ENDIF CALL ROT2(VEC,THETA0-THETA,RM2) CALL MULT3(RM2,RM1,RM3) C**** Now rotate the target through angle BETA-BETA0 around a1 C (i.e., rotate the lab through angle BETA0-BETA around a1) CALL ROT2(A1,BETA0-BETA,RM1) CALL MULT3(RM1,RM3,RM2) C C**** RM2 is now the full rotation matrix C C**** Now rotate all required vectors: C Incident propagation vector: VEC(1)=1. VEC(2)=0. VEC(3)=0. CALL PROD3(RM2,VEC,EN0R) C Incident polarization vectors: CALL PROD3C(RM2,CXE01,CXE01R) CALL PROD3C(RM2,CXE02,CXE02R) C Scattering vectors: CALL PROD3V(RM2,ENSC,AKSR,NSCAT,MXSCA) DO 50 I=1,NSCAT AKSR(1,I)=AK1*AKSR(1,I) AKSR(2,I)=AK1*AKSR(2,I) AKSR(3,I)=AK1*AKSR(3,I) 50 CONTINUE C Scattering polarization vectors: CALL PROD3V(RM2,EM1,EM1R,NSCAT,MXSCA) CALL PROD3V(RM2,EM2,EM2R,NSCAT,MXSCA) RETURN END SUBROUTINE SCAT(AK,AKS,EM1,EM2,E02,CMDTRQ, & CBKSCA,CSCA,CSCAG,CTRQIN,CTRQSC, & CXE,CXE00,CXF1L,CXF2L,CXP,CXSCR1,CXSCR2,CXSCR3,CXSCR4, & ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR, & PHIN,SCRRS1,SCRRS2,THETAN,IXYZ) C C Arguments: C CHARACTER CMDTRQ*6 REAL CBKSCA,CSCA,E02 INTEGER ICTHM,IPHIM,MXN3,MXNAT,MXSCA,NAT,NAT3,NDIR INTEGER*2 IXYZ(MXNAT,3) COMPLEX CXE(NAT,3),CXE00(3),CXF1L(MXSCA),CXF2L(MXSCA),CXP(NAT,3), & CXSCR1(MXN3),CXSCR2(MXN3),CXSCR3(MXN3),CXSCR4(MXNAT,3) REAL AK(3),AKS(3,MXSCA),CSCAG(3),CTRQIN(3),CTRQSC(3), & EM1(3,MXSCA),EM2(3,MXSCA), & PHIN(MXSCA),SCRRS1(NAT,3),SCRRS2(MXNAT),THETAN(MXSCA) C C Local variables: C COMPLEX CXI,CXSCL1,CXSCL2,CXSCL3 REAL AK2,AK3,AKK,CONST,COSPHI,COSTH,DCOSTH,DPHI,EINC,ES2, & PI,PHI,RRR,SINPHI,SINTH,TERM,W,WTH INTEGER ICOSTH,IPHI,IPHM,J,K,ND COMPLEX CXBS(3),CXES(3),CXTRM(3) REAL AKS0(3),AKS1(3),AKS2(3),AKSN(3),CTRQSCTF(3),XYZCM(3) C C Intrinsic functions: C INTRINSIC COS,EXP,MIN,SIN,SQRT,CONJG,REAL C C Data statements: C DATA CXI/(0.,1.)/ C*********************************************************************** C C SCAT computes energy radiated by array of oscillating dipoles C and corresponding scattering properties if oscillations are in C response to incident E field. C C Given: C AK(1-3) = KX,KY,KZ = components of incident k-vector in TF C AKS(3,MXSCA)=scattered k vectors in TF C CMDTRQ = 'DOTORQ' to calculate torques C = 'NOTORQ' to skip calculation of torques C EM1(3,MXSCA)=pol.vector 1 for each scattering direction in TF C EM2(3,MXSCA)=pol.vector 2 for each scattering direction in TF C E02 = |E_0|^2 , where E_0 = incident complex E field C NAT = number of dipoles C NAT3 = 3*NAT C IXYZ(NAT,1-3) = x/d, y/d, z/d (integers) for dipoles 1-NAT C CXE(1-NAT,1-3) = components of E field at each dipole at t=0 C CXE00(1-3) = (complex) reference polarization vector C CXP(1-NAT,1-3) = components of polarization of each dipole C ICTHM = number of THETA values for calculation of CSCA and G C (ICTHM must be odd; recommend ICTHM > 5*(2*pi*a/lambda) C for accuracy) C IPHIM = number of PHI values for calculation of CSCA and G C NDIR = number of directions at which scattering matrix elements C are to be computed C THETAN(1-NDIR) : scattering angle (radians); C cos(THETAN) = n0 dot n C *** BTD 98.11.18: THETAN is not used *** C PHIN(1-NDIR) : scattering angle (radians); C *** BTD 98.11.18: PHIN is not used *** C C Returns: C CBKSCA=differential scattering cross section for theta=pi C (lattice units) C CSCA = scattering cross section (lattice units) C CSCAG(1) = CSCA* , where theta=scattering angle C CSCAG(2) = CSCA* C CSCAG(3) = CSCA* C CTRQIN(1-3)=cross section for torque on grain due to incident C E and B fields, C for torque in directions x,y,z assuming incident C radiation is along x direction, C and "reference pol state" is in y direction. C CTRQSC(1-3)=cross section for torque on grain due to scattered C E and B fields, C for torque in directions x,y,z assuming incident C radiation is along x direction, C and "reference pol state" is in y direction C Here "torque cross section" is defined so that C torque(1-3)=(momentum flux)*(1/k)*C_torque(1-3) C Total torque cross section=CTRQIN+CTRQSC C CXF1L(NDIR)=f_{1L} for direction NDIR C CXF2L(NDIR)=f_{2L} for direction NDIR C where f_{1L} connects incident polarization L to C outgoing theta polarization (E in scattering plane) C f_{2L} connects incident polarization L to outgoing phi C polarization (E perpendicular to scattering plane) C This is same f_{ml} notation as used by Draine (1988; C Ap.J., 333, 848). C C Notes: C AKS0(1-3)=k_s=scattered k vector (lattice units) C AKSN(1-3)=k_s/|k_s|=unit vector in scattering direction C Scratch vectors introduced to permit vectorization: C SCRRS1 = real vector of length.GE.3*NAT C SCRRS2 = real vector of length.GE.NAT C CXSCR1 = complex scratch array of length 3*MXNAT C CXSCR2 = complex scratch array of length 3*MXNAT C CXSCR3 = complex scratch array of length 3*MXNAT C CXSCR4 = complex scratch array of length 3*MXNAT C C B.T.Draine, Princeton Univ. Obs., 87.01.04 C History: C 88.06.29 (BTD): modified C 89.11.28 (BTD): modified C 90.11.02 (BTD): modified to enable use with vacuum sites. C 90.11.20 (BTD): corrected error in calculation of f_ml (missing C factor of k). C 90.12.03 (BTD): changed ordering of XYZ0 and CXP C eliminated cxe0 from argument list C 90.12.10 (BTD): replaced XYZ0 with IXYZ0 to reduce memory use C 91.08.14 (BTD): add CBKSCA to argument list for SCAT C add code to eliminate sum over phi when theta=0 or pi C add code to calculate CBKSCA C 94.04.04 (BTD): removed superfluous CMPLX( ) operation from 2 lines C in DO 200 loop. Also added new variable TERM to move C some computation out of last DO loop. C 95.06.14 (BTD): modified to compute moments C and of scattered radiation C added code to compute torque on grain relative to C grain centroid XYZCM(1-3) C NOTE: handling of grain centroid is not as efficient C as it should be. Can be recoded to reduce number C of computations. C 95.06.15 (BTD): modified to compute absorption contribution to C torque on grain. C added CXE and CTRQIN to argument list C 95.06.15 (BTD): corrected errors found by Joseph Weingartner C 95.06.16 (BTD): more corrections C 95.06.20 (BTD): added argument CMDTRQ to allow torque calculations C to be skipped if not desired. C 95.06.22 (BTD): added factor CXI in terms appearing in evaluation C of scattered torque C 95.06.29 (BTD): Added loop to compute AKS0(1-3) for special cases C theta-0 and 180. Error found by J. Weingartner. C 95.07.20 (BTD): Add scratch vector SCRRS2 C Add scratch array CXSCR4 to argument list C Rewrite to make torque computation more efficient: C move as many calculations as possible out of loops C over scattering directions: C SCRRS1 now used to store coords relative to centroid C SCRRS2 now used to store (k_s dot r_j) C CXSCR4 is needed to store (p_j cross r_j) C 95.07.21 (BTD): Rewrite (again!) to make the correspondence C to notation of Draine & Weingartner more apparent. C This involved changing various definitions by C factors of k, and change of sign of CXBS C 95.07.24 (BTD): Corrected errors in evaluation of both scattering C torque and incident torque C 95.07.28 (BTD): Reenabled computation of CBKSCA C 98.11.18 (BTD): Edited comments C C Copyright (C) 1993,1994,1995,1998 B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** PI=4.*ATAN(1.) C C ICTHM,IPHIM determine angular grid used for calculation of total C scattering cross section and G= C ICTHM=number of theta values for which scattering is evaluated C (0 to pi) C IPHIM=number of phi values for which scattering is evaluated C (0 to 2*pi) C Note: ICTHM must be odd as Simpson's rule integration is employed. C Check whether ICTHM is odd, and increase by 1 if not. C IF(ICTHM.EQ.2*(ICTHM/2))ICTHM=ICTHM+1 C C Incident k vector AK defines one axis (for spherical integration) C Use (real part of) reference polarization vector CXE00 to define C second coordinate axis. C Construct third coordinate axis by taking cross product of first C two axes. C AKS1, AKS2 = vectors of length equal to AK lying along second and C third coordinate axes. C AK2=AK(1)*AK(1)+AK(2)*AK(2)+AK(3)*AK(3) AKK=SQRT(AK2) AK3=AKK*AK2 AKS1(1)=REAL(CXE00(1)) AKS1(2)=REAL(CXE00(2)) AKS1(3)=REAL(CXE00(3)) C C EINC = magnitude of (real part of) reference electric field C EINC=SQRT(AKS1(1)**2+AKS1(2)**2+AKS1(3)**2) IF(EINC.LE.0.)CALL ERRMSG('FATAL','SCAT', & 'CXE00 IS PURE IMAGINARY!') RRR=AKK/EINC AKS1(1)=RRR*AKS1(1) AKS1(2)=RRR*AKS1(2) AKS1(3)=RRR*AKS1(3) AKS2(1)=(AK(2)*AKS1(3)-AK(3)*AKS1(2))/AKK AKS2(2)=(AK(3)*AKS1(1)-AK(1)*AKS1(3))/AKK AKS2(3)=(AK(1)*AKS1(2)-AK(2)*AKS1(1))/AKK C DCOSTH=2./REAL(ICTHM-1) DPHI=2.*PI/REAL(IPHIM) CSCA=0. DO K=1,3 CSCAG(K)=0. CTRQSCTF(K)=0. ENDDO C C CONST = multiplicative factor entering weights C factor 2/(ICTHM-1) = d(cos(theta)) C factor 2*pi/IPHIM = d(phi) C factor 3 in denom from Simpson's rule (later mult. by 1,4,2..) C Use Simpson's rule integration over d(cos(theta)) C CONST=4.*PI/ (3.*REAL((ICTHM-1)*IPHIM)) C C Note that for cos(theta)=-1 or 1 we do not need to sum over phi, C so treat these two directions separately C IF(CMDTRQ.EQ.'DOTORQ')THEN C C*********************************************************************** C Additional quantities required for torque computations C XYZCM(1-3) = coordinates of grain centroid C SCRRS1(J,K)= r_j = coordinates of dipole J relative to centroid. C CXSCR3(J) = r_j dot p_j C CXSCR4(J,1-3) = p_j cross r_j C DO K=1,3 XYZCM(K)=0. DO J=1,NAT XYZCM(K)=XYZCM(K)+IXYZ(J,K) ENDDO XYZCM(K)=XYZCM(K)/NAT ENDDO DO K=1,3 DO J=1,NAT SCRRS1(J,K)=IXYZ(J,K)-XYZCM(K) ENDDO ENDDO DO J=1,NAT CXSCR3(J)=(0.,0.) DO K=1,3 CXSCR3(J)=CXSCR3(J)+SCRRS1(J,K)*CXP(J,K) ENDDO ENDDO DO J=1,NAT CXSCR4(J,1)=CXP(J,2)*SCRRS1(J,3)-CXP(J,3)*SCRRS1(J,2) CXSCR4(J,2)=CXP(J,3)*SCRRS1(J,1)-CXP(J,1)*SCRRS1(J,3) CXSCR4(J,3)=CXP(J,1)*SCRRS1(J,2)-CXP(J,2)*SCRRS1(J,1) ENDDO C*********************************************************************** ELSE DO K=1,3 DO J=1,NAT SCRRS1(J,K)=IXYZ(J,K) ENDDO ENDDO ENDIF C C Proceed to sum over scattering angles: C DO ICOSTH=1,ICTHM IPHM=IPHIM WTH=4. IF(ICOSTH.GT.2*(ICOSTH/2))WTH=2. C C In special case ICOSTH=1 or ICTHM, we have theta=0 or pi, C and the angle phi becomes irrelevant. In this case we C need only one value of phi, so we set IPHM=1 and adjust weights C accordingly. C IF(ICOSTH.EQ.1.OR.ICOSTH.EQ.ICTHM)THEN WTH=IPHIM IPHM=1 ENDIF COSTH=1.-REAL(ICOSTH-1)*DCOSTH SINTH=SQRT(1.-MIN(COSTH*COSTH,1.)) C C Evaluate total weight factor W C W=CONST*WTH DO IPHI=1,IPHM PHI=DPHI*(REAL(IPHI)-0.5) SINPHI=SIN(PHI) COSPHI=COS(PHI) C C Compute scattered k vector AKS0(1-3) C Initialize CXES(1-3) = scattered electric field C DO K=1,3 AKS0(K)=COSTH*AK(K)+SINTH* & (COSPHI*AKS1(K)+SINPHI*AKS2(K)) CXES(K)=(0.,0.) ENDDO C C Compute normalized scattering vector = nhat C DO K=1,3 AKSN(K)=AKS0(K)/AKK ENDDO C C Now compute CXES(1-3) = sum_j[p_j-nhat(nhat dot p_j)]exp[-ik_s dot x_j] C =(r/k*k)*exp(-ikr+iwt)*E(k_s,r,t) C where E(k_s,r,t) is complex electric field vector propagating C with wave vector k_s, at distance r from origin C Notation: C If CMDTRQ='DOTORQ', then r_j = position relative to centroid XYZCM C If CMDTRQ='NOTORQ', then r_j = original lattice coords, and XYZCM=0. C Note that since we are not concerned with overall phase shifts for C evaluation of forces and torques, we do NOT need to worry about change C in phase resulting from using two different definitions of r_j when C computing phase factors exp[-i k_s dot r_j] C DO J=1,NAT C C CXSCR1(J) = nhat dot p(j) C CXSCR1(J)=(0.,0.) DO K=1,3 CXSCR1(J)=CXSCR1(J)+AKSN(K)*CXP(J,K) ENDDO C C SCRRS2(J) = nhat dot r(j) C SCRRS2(J)=0. DO K=1,3 SCRRS2(J)=SCRRS2(J)+SCRRS1(J,K)*AKSN(K) ENDDO C C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j C CXSCR2(J)=EXP(-CXI*AKK*SCRRS2(J)) C C CXES(1-3) = sum_j[p_j-nhat(nhat dot p_j)]exp[-ik_s dot x_j] C DO K=1,3 CXES(K)=CXES(K)+(CXP(J,K)-AKSN(K)*CXSCR1(J))* & CXSCR2(J) ENDDO ENDDO C IF(CMDTRQ.EQ.'DOTORQ')THEN C*********************************************************************** C Additional terms required to compute torque on grain C Have already computed C CXES(1-3)=sum_j[p_j - k_s*(k_s dot p_j)/|k|^2]exp[-ik_s dot x_j] C where p_j = polarization of dipole j C x_j = position of dipole j rel. to CM C k_s = scattered k vector C C CXBS = -nhat cross CXES = -(r/k*k)*exp(-ikr+iwt)*B(k_s,r,t) C CXBS(1)=AKSN(3)*CXES(2)-AKSN(2)*CXES(3) CXBS(2)=AKSN(1)*CXES(3)-AKSN(3)*CXES(1) CXBS(3)=AKSN(2)*CXES(1)-AKSN(1)*CXES(2) C C CXSCL1 = sum_j p_j dot (r_j cross nhat)*exp(-i*k_s dot r_j) C = nhat dot sum_j (p_j cross r_j)*exp(-i*k_s dot r_j) C recall that have previously evaluated C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j C CXSCR4(J,K)=p_j cross r_j C CXSCL1=(0.,0.) DO K=1,3 CXTRM(K)=(0.,0.) ENDDO DO K=1,3 DO J=1,NAT CXTRM(K)=CXTRM(K)+CXSCR2(J)*CXSCR4(J,K) ENDDO CXSCL1=CXSCL1+AKSN(K)*CXTRM(K) ENDDO C C CXSCL2 = sum_j[(r_j dot p_j) - C (nhat dot p_j)*((nhat dot r_j) + 2i/k_s)]* C exp(-i*k_s dot r_j) C Have previously computed C CXSCR1(J) = nhat dot p(j) C CXSCR2(J) = exp(-i*k_s dot r_j) factor for atom j C CXSCR3(J) = r_j dot p_j C SCRRS2(J) = nhat dot r(j) C CXSCL2=(0.,0.) DO J=1,NAT CXSCL2=CXSCL2+CXSCR2(J)*( & CXSCR3(J)-CXSCR1(J)*(SCRRS2(J)+2.*CXI/AKK)) ENDDO C C Note following notational correspondence: C present Draine & Weingartner (1996) C ------- --------------------------- C CXES = V_E C CXBS = V_B C CXSCL1 = S_B C CXSCL2 = S_E C C*********************************************************************** ENDIF C C Have completed computation of vector CXES(1-3) and terms to calculate C angular momentum flux for scattering direction (THETA,PHI). C CTRQSCTF(K) is "torque cross section" due to "scattered radiation" C resulting in torque in lattice direction k (i.e., in "target C frame". After CTRQSCTF has been evaluated, we will compute CTRQSC in C the "scattering frame", where the radiation is incident along x, C and the "reference polarization state" is along the y direction. C ES2=0.0 DO K=1,3 ES2=ES2+CXES(K)*CONJG(CXES(K)) ENDDO CSCA=CSCA+W*ES2 CSCAG(1)=CSCAG(1)+W*COSTH*ES2 CSCAG(2)=CSCAG(2)+W*SINTH*COSPHI*ES2 CSCAG(3)=CSCAG(3)+W*SINTH*SINPHI*ES2 IF(CMDTRQ.EQ.'DOTORQ')THEN DO K=1,3 CTRQSCTF(K)=CTRQSCTF(K)- & W*REAL(CONJG(CXBS(K))*CXSCL2+CONJG(CXES(K))*CXSCL1) ENDDO ENDIF ENDDO C C Save backscattering cross section: C IF(ICOSTH.EQ.ICTHM)CBKSCA=AK2*AK2*ES2/E02 ENDDO C C*********************************************************************** C CSCA=AK2*AK2*CSCA/E02 C C Convert CSCAG and CTRQSCTF to cross sections measured in lattice units: C DO K=1,3 CSCAG(K)=AK2*AK2*CSCAG(K)/E02 ENDDO IF(CMDTRQ.EQ.'DOTORQ')THEN C C*********************************************************************** C Compute the two contributions to the torque on the target: C (1) torque due to E and B fields of incident radiation (CTRQIN) C (2) torque due to E and B fields of scattered radiation (CTRQSC) C These are first computed in the target frame, then transformed to C the lab frame (frame where incident radiation propagates along C the x axis, and "reference pol. state 1" is along the y axis). C DO K=1,3 CTRQSCTF(K)=AK2*AK2*CTRQSCTF(K)/E02 ENDDO C C CTRQSCTF(1-3) is the vector torque cross section/|k| measured in the C target frame. We now compute CTRQSC(1-3)=the vector torque cross C section in the "scattering frame", where the radiation is incident along C the x axis, and the "reference polarization state" (at t=0) is C linearly polarized along the y axis. C Note that the following computation introduces an additional factor C of |k|. C CTRQSC(1)=AK(1)*CTRQSCTF(1)+AK(2)*CTRQSCTF(2)+AK(3)*CTRQSCTF(3) CTRQSC(2)=AKS1(1)*CTRQSCTF(1)+AKS1(2)*CTRQSCTF(2)+ & AKS1(3)*CTRQSCTF(3) CTRQSC(3)=AKS2(1)*CTRQSCTF(1)+AKS2(2)*CTRQSCTF(2)+ & AKS2(3)*CTRQSCTF(3) C C Now compute vector torque cross section due to incident radiation: C DO K=1,3 CTRQIN(K)=0. ENDDO C C Redefine scratch vector CXSCR1: C CXSCR1(J)=conjg(p_j) dot E_0j at t=0 C DO J=1,NAT CXSCR1(J)=0. DO K=1,3 CXSCR1(J)=CXSCR1(J)+CONJG(CXP(J,K))*CXE(J,K) ENDDO ENDDO CXSCL1=0. CXSCL2=0. CXSCL3=0. C C Compute CXSCL1,CXSCL2,CXSCL3=components of the torque due to C incident E and B fields in the target frame. C 95.06.22 (BTD): added factor CXI in r cross k term appearing in C evaluations of CXSCL1,CXSCL2,CXSCL2 C 95.07.21 (BTD): corrected sign mistake (+CXI -> -CXI) C 95.07.24 (BTD): changed def. of CXSCR1 to agree with paper, C changed calc. of CXSCLn to more closely mirror paper, C and changed AKS0 to AK C DO J=1,NAT CXSCL1=CXSCL1+CONJG(CXP(J,2))*CXE(J,3)- & CONJG(CXP(J,3))*CXE(J,2)-CXI* & (AK(2)*SCRRS1(J,3)-AK(3)*SCRRS1(J,2))*CXSCR1(J) CXSCL2=CXSCL2+CONJG(CXP(J,3))*CXE(J,1)- & CONJG(CXP(J,1))*CXE(J,3)-CXI* & (AK(3)*SCRRS1(J,1)-AK(1)*SCRRS1(J,3))*CXSCR1(J) CXSCL3=CXSCL3+CONJG(CXP(J,1))*CXE(J,2)- & CONJG(CXP(J,2))*CXE(J,1)-CXI* & (AK(1)*SCRRS1(J,2)-AK(2)*SCRRS1(J,1))*CXSCR1(J) ENDDO C C Compute CTRQIN(1-3)=2*components of torque due to incident E and B C fields (in lab frame), expressed as a cross section: C CTRQIN(1)=AK(1)*REAL(CXSCL1)+AK(2)*REAL(CXSCL2)+ & AK(3)*REAL(CXSCL3) CTRQIN(2)=AKS1(1)*REAL(CXSCL1)+AKS1(2)*REAL(CXSCL2)+ & AKS1(3)*REAL(CXSCL3) CTRQIN(3)=AKS2(1)*REAL(CXSCL1)+AKS2(2)*REAL(CXSCL2)+ & AKS2(3)*REAL(CXSCL3) DO K=1,3 CTRQIN(K)=4.*PI*CTRQIN(K)/E02 ENDDO ENDIF C C*********************************************************************** C IF(NDIR.LE.0)RETURN C C*** Calculate the 2 matrix elements FM1 for selected directions C F1L(NDIR) is to the polarization state e1 (E in scattering plane) C F2L(NDIR) is to the polarization state e2 (E perp. to scatt. plane) C C IMPORTANT NOTE: If the incident radiation is not linearly polarized, C with zero imaginary component to the polarization vector, then the C F1L and F2L computed here are NOT the "usual" f_ml, which C describe scattering from linearly polarized states l to linearly C polarized states m. C If one wishes to reconstruct the "usual" f_ml even when computing C with elliptically polarized light, then one must add additional C code to the routine SCAT to reconstruct the f_ml from the F1L and C F2L computed here. C TERM=AK3/SQRT(E02) DO ND=1,NDIR CXF1L(ND)=(0.0,0.0) CXF2L(ND)=(0.0,0.0) DO J=1,NAT CXSCR2(J)=(0.0,0.0) CXSCR3(J)=(0.0,0.0) DO K=1,3 CXSCR2(J)=CXSCR2(J)+CXP(J,K)*EM1(K,ND) CXSCR3(J)=CXSCR3(J)+CXP(J,K)*EM2(K,ND) ENDDO ENDDO DO J=1,NAT CXSCR1(J)=EXP(-CXI*(AKS(1,ND)*IXYZ(J,1)+ & AKS(2,ND)*IXYZ(J,2)+ & AKS(3,ND)*IXYZ(J,3))) CXF1L(ND)=CXF1L(ND)+CXSCR2(J)*CXSCR1(J) CXF2L(ND)=CXF2L(ND)+CXSCR3(J)*CXSCR1(J) ENDDO C C Units: Dipole polarizations are in units of [E]*d**3, where d=lattice C spacing, and [E]=dimensions of E field C Therefore, up to this point CXF1L, CXF2L are in units C of [E]*d**3 C Now multiply by AK3/SQRT(E02) (units of 1/([E]d**3)) to get C dimensionless result. C CXF1L(ND)=CXF1L(ND)*TERM CXF2L(ND)=CXF2L(ND)*TERM ENDDO RETURN C END SUBROUTINE SCAVEC(MXSCA,NSCAT,THETAN,PHIN,CXE01,ENSC,EM1,EM2) C*********************************************************************** C Given: C MXSCA=dimension information for ENSC,PHIN,THETAN C NSCAT=number of scattering directions C THETAN(1-NSCAT)=scattering angles theta C PHIN(1-NSCAT)=scattering angles phi C CXE01(1-3)=complex polarization vector 1 (phi=0 direction C is defined by x,y plane (Lab Frame), where incident C radiation propagates along the x axis. C Returns: C ENSC(1-3,1-NSCAT)=scattering vectors in Lab Frame C EM1(1-3,1-NSCAT)=scattered pol vectors parallel to scat. plane C in Lab Frame C EM2(1-3,1-NSCAT)=scattered pol vectors perp. to scat. plane C in Lab Frame C It is assumed that incident propagation vector is in x-direction C in Lab Frame C C History: C 96.11.06 (BTD): Changed definition of scattering angle phi C Previously, phi was measured from plane containing C incident k vector (i.e., Lab x-axis) and Re(CXE01) C Henceforth, phi is measured from Lab x,y plane. C C Copyright (C) 1993,1996, B.T. Draine and P.J. Flatau C This code is covered by the GNU General Public License. C*********************************************************************** C C Arguments: INTEGER MXSCA,NSCAT COMPLEX CXE01(3) REAL EM1(3,MXSCA), EM2(3,MXSCA), ENSC(3,MXSCA), PHIN(MXSCA), & THETAN(MXSCA) C C Local variables: INTEGER I REAL COSPHI,COSTHE,SINPHI,SINTHE C C*********************************************************************** DO 1000 I=1,NSCAT COSPHI=COS(PHIN(I)) SINPHI=SIN(PHIN(I)) COSTHE=COS(THETAN(I)) SINTHE=SIN(THETAN(I)) ENSC(1,I)=COSTHE ENSC(2,I)=SINTHE*COSPHI ENSC(3,I)=SINTHE*SINPHI EM1(1,I)=-SINTHE EM1(2,I)=COSTHE*COSPHI EM1(3,I)=COSTHE*SINPHI EM2(1,I)=0. EM2(2,I)=-SINPHI EM2(3,I)=COSPHI 1000 CONTINUE RETURN END