C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c KRA - substitution coordinates from single isotopic substitution in c a linear/symmetric/asymmetric top from Kraitchman's equations c C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ c c - The program uses single isotopic equations straight out of c Gordy&Cook c - Experimental errors in rotational constants are propagated into c errors in coordinates c - The much more realistic errors according to the Costain criterion c are then added to the errors above, c - For planar molecules suitable planarity conditions can be built in c to use any two out of the three measured constants c c C Version 14d.X.2002 ----- Zbigniew KISIEL ----- C C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C | http://info.ifpan.edu.pl/~kisiel/prospe.htm | C_________________________/-------------------------------------------------- c c Modification history: c c ca 95: Creation c 23.10.97: Addition of reference data and modifications to output c 30.05.99: Debugging of input c 1.04.01: Merging with symmetric top equations from another version c 14.10.02: Some output mods C_____________________________________________________________________________ c c Input data set is to contain: c c line0: comment line c line1: i j c line1a: k l c line2: A, B, C \ c line3: dA, dB, dC |- for the parent species c line4: Mass / c line5: A, B, C \ c line6: dA, dB, dC |- for a substituted species c line7: dMass / c . c . c c - Line1 selects the type of top: c c i=0: defines that the parent is a diatomic, linear or a symmetric top c (off-axis asymmetric substitution is allowed for a symmetric top) c If i=0 then the value of j is irrelevant but is still necessary on c input so it is best to also set it to zero. c c i>0: defines that the parent is an asymmetric top, in which case i and j c can take values out of 1,2,3 (i.e. a,b,c). i cannot be equal to j, c and (i,j) selects the inertial plane for planar coordinate c calculation. c One of i and j can be negative in which case the respective moment c of inertia is calculated using the planar relation defined in line1a c c - Line1a is only necessary if one of i,j is negative and defines the c planar relation for calculation of that moment, for example c if i,j = -1 2 c then k,l = 3 -2 define I.a = I.c - I.b c c - Line 2: rotational constants for the parent species (MHz). If the c species is not an asymmetric top then set A to a value large enough to c pretend that it is infinity, say to 9999999.0 c c - Line 3: errors in rotational constants for the parent (MHz) c c - Line 4: mass of the parent isotopic species in u c c - Lines 5,6: rotational constants and errors for a substituted species c c - Line 7: the mass difference between the mass of substituted species and c that of the parent species (u), see the table below for some common values. c c - Lines 5-7 can be repeated as many times as necessary, c c NOTE: c c - Lines that begin with a recognised comment flag ie. c with either one of c,C,#,!,% are ignored so that c comment lines can be inserted anywhere in the data file with c the exception of the first line, which is always treated as c comment describing the data set c c - The X,Y,Z coordinates in the output are consistent with the c x,y,z indices in the equations in Gordy&Cook c c - Remember that all substitution coordinates determined with Kraitchman's c equations arise from a square root operation so that their sign is c indeterminate. For convenience KRA prints all coordinates positive, and c not with +- sign in front of each. c It is, of course, necessary to assign signs to coordinates for their c use in determination of internal structural parameters, as for example c with EVAL. c C_____________________________________________________________________________ C C References: C C [1] W.Gordy and R.L.Cook, "Microwave Molecular Spectra", Wiley, C New York, 1984. c The relevant relations are: c linear: Eq.13.46 c off axis symm.: Eq.13.56-13.57 c planar asymm: Eq.13.64-13.65 c general asymm. Eq.13.71-13.73 c C [2] C.C.Costain, Trans.Am.Crystallogr.Assoc. 2,157-164(1966). C_____________________________________________________________________________ C C C M.parent (u) M.isot (u) Delta.M (u) C c 1.H = 1.0078252 2.H = 2.0141022 1.0062770 C 12.C = 12. 13.C = 13.0033544 1.0033544 c 14.N = 14.0030744 15.N = 15.0001077 0.9970333 c 16.O = 15.9949150 18.O = 17.999160 2.004245 c 17.O = 16.999133 1.004218 c 32.S = 31.9720737 33.S = 32.971462 0.9993883 c 34.S = 33.967867 1.9957933 c 35.Cl = 34.968851 37.Cl = 36.965899 1.997048 c 79.Br = 78.918329 81.Br = 80.916292 1.997963 c C C_____________________________________________________________________________ C c IMPLICIT REAL*8(A-H,O-Z) COMMON RIX,RIY,RIZ,RIX1,RIY1,RIZ1,RM, * erix,eriy,eriz,erix1,eriy1,eriz1 common /moment/parent(3,2),isot(3,2),iflagx,iflagy,ifirst,isecnd PARAMETER (CONV=505379.005D0) integer iflagx,iflagy,ifirst,isecnd real*8 isot character*2 ixx,iyy,izz,ixxp,iyyp,irecal character filnam*35,type*13,coment*78,fsign*1,ssign*1,errmes*50 CHARACTER*3 AX(3) c equivalence * (parent(1,1),ria),(parent(2,1),rib),(parent(3,1),ric), * (parent(1,2),eria),(parent(2,2),erib),(parent(3,2),eric), * (isot(1,1),ria1),(isot(2,1),rib1),(isot(3,1),ric1), * (isot(1,2),eria1),(isot(2,2),erib1),(isot(3,2),eric1) DATA AX/'a','b','c'/ isym=0 iflagx=2 iflagy=3 C 2 CONTINUE write(*,'(1x//)') WRITE(*,3344) 3344 FORMAT(' ',76(1H_)/' |',T79,'|'/ * ' | KRA - SINGLE ISOTOPIC SUBSTITUTION - Various ', * 'permutations',T79,'|'/ * ' | of Kraitchman''s equat', * 'ions for symmetric/asymmetric tops',T79,'|'/ * ' |',76(1H_),'|'/' version 14d.X.2002',T64,'Zbigniew KISIEL'/) C 120 write(*,'(1x/'' Name of data file: '',$)') read(*,'(a)',err=120)filnam open(4,file=filnam,status='old',err=120) read(4,'(a)',err=120)coment open(3,file='kra.out',status='unknown') c write(3,3344) write(3,'(1x,78(1h-)/1x,a/1x,78(1h-)/)')coment c c c...data for parent species c 1121 errmes='Problems reading i,j (internally IFLAGX,IFLAGY)' 2002 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 2002 read(4,*,err=2000,end=99)iflagx,iflagy if(iflagx.eq.0)then isym=1 goto 121 endif write(errmes,'(2i7)')iflagx,iflagy errmes='Illegal i,j = '//errmes(1:14) if(iflagx.lt.0.and.iflagy.lt.0)goto 2000 if(iabs(iflagx).lt.1.or.iabs(iflagx).gt.3)goto 2000 if(iabs(iflagy).lt.1.or.iabs(iflagy).gt.3)goto 2000 if(iflagx.eq.iflagy)goto 2000 c if(iflagx.lt.0.or.iflagy.lt.0)then errmes='Problems reading k,l (internally IFIRST,ISECND)' 2003 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 2003 read(4,*,err=2000,end=99)ifirst,isecnd write(errmes,'(2i7)')ifirst,isecnd errmes='Illegal k,l = '//errmes(1:14) if(iabs(ifirst).lt.1.or.iabs(ifirst).gt.3)goto 2000 if(iabs(isecnd).lt.1.or.iabs(isecnd).gt.3)goto 2000 if(iabs(ifirst).eq.iabs(isecnd))goto 2000 fsign=' ' ssign='+' if(ifirst.lt.0)fsign='-' if(isecnd.lt.0)ssign='-' if(iflagx.lt.0)write(3,37)ax(iabs(iflagx)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)), * ax(iabs(iflagy)) if(iflagy.lt.0)write(3,38)ax(iabs(iflagx)), * ax(iabs(iflagy)), * fsign,ax(iabs(ifirst)),ssign,ax(iabs(isecnd)) else write(*,*)iflagx,iflagy write(3,36)ax(iabs(iflagx)),ax(iabs(iflagy)) endif c 36 format(1x/' Planar calculation will be made from I.',a1, * ' and I.',a1/) 37 format(1x/' Planar calculation will be made from I.',a1, * ' = ',a1,'I.',a1,' ',a1,' I.',a1,' and I.',a1/) 38 format(1x/' Planar calculation will be made from I.',a1, * ' and I.',a1,' = ',a1,'I.',a1,' ',a1,' I.',a1/) c errmes='Problem reading A,B,C' 121 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 121 read(4,*,err=2000,end=99)a,b,c c errmes='Problem reading dA,dB,dC' 122 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 122 read(4,*,err=2000,end=99)era,erb,erc c errmes='Problem reading mass of parent' 123 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 123 read(4,*,err=2000,end=99)tm c RIA=CONV/A RIB=CONV/B RIC=CONV/C eria=conv*era/a**2 erib=conv*erb/b**2 eric=conv*erc/c**2 defect=RIC-RIA-RIB rkappa=(2.d0*b-a-c)/(a-c) APLUSB=(A+B)/C PA=0.5D0*(RIB+RIC-RIA) PB=0.5D0*(RIA+RIC-RIB) PC=0.5D0*(RIA+RIB-RIC) write(3,'(1x,78(1h-))') C c...data for isotopic species (looped over as many times as in file) c errmes='Problem reading isotopic A,B,C' 124 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 124 READ(4,*,err=2000,end=99)A1,B1,C1 if(a1.le.0.0)goto 2000 c errmes='Problem reading isotopic dA,dB,dC' 125 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 125 read(4,*,err=2000,end=99)era1,erb1,erc1 if(era1.le.0.0)goto 2000 c errmes='Problem reading isotopic mass change' 126 call test(4,ierr) <--- TEST if(ierr.eq.1)goto 126 READ(4,*,err=2000,end=99)DM if(dm.eq.0.0.or.dabs(dm).gt.100.)goto 2000 RM=TM*DM/(TM+DM) C RIA1=CONV/A1 RIB1=CONV/B1 RIC1=CONV/C1 eria1=conv*era1/a1**2 erib1=conv*erb1/b1**2 eric1=conv*erc1/c1**2 defec1=RIC1-RIA1-RIB1 rkapp1=(2.d0*b1-a1-c1)/(a1-c1) APLUS1=(A1+B1)/C1 PA1=0.5D0*(RIB1+RIC1-RIA1) PB1=0.5D0*(RIA1+RIC1-RIB1) PC1=0.5D0*(RIA1+RIB1-RIC1) C c...trace of data that has been input c write(3,'(1x,78(1h-))') WRITE(*,10) 10 FORMAT(1X/' ======> The common species:'/) WRITE(*,33)' A, B, C = ',A,B,C write(*,33)' dA, dB, dC = ',era,erb,erc WRITE(*,33)' IA, IB, IC = ',RIA,RIB,RIC WRITE(*,33)'dIA,dIB,dIC = ',erIA,eRIB,eRIC c WRITE(*,33)'PA,PB,PC = ',PA,PB,PC c WRITE(*,33)'In.Defect= ',defect c WRITE(*,33)'Kappa = ',rkappa c WRITE(*,33)'A+B /C = ',APLUSB WRITE(*,35)'Mass = ',TM WRITE(*,15) 15 FORMAT(1X/' ======> The isotopic species:'/) WRITE(*,33)' A, B, C = ',A1,B1,C1 write(*,33)' dA, dB, dC = ',era1,erb1,erc1 WRITE(*,33)' IA, IB, IC = ',RIA1,RIB1,RIC1 WRITE(*,33)'dIA,dIB,dIC = ',eRIA1,eRIB1,eRIC1 c WRITE(*,33)'PA,PB,PC = ',PA1,PB1,PC1 c WRITE(*,33)'In.Defect= ',defec1 c WRITE(*,33)'Kappa = ',rkapp1 c WRITE(*,33)'A+B /C = ',APLUS1 WRITE(*,35)'Mass change = ',dm write(*,35)'Total mass = ',TM+DM 33 FORMAT(1X,A,3F20.8) 35 format(1x,a,f14.8) c c c...........KRAITCHMAN calculation............... c WRITE(*,43) 43 FORMAT(/' KRAITCHMAN''S RESULTS '/) c c...Linear-symmetric top calculation c if(isym.eq.1)then if(rib.eq.ric)then xx=c yy=b zz=a erxx=erc eryy=erb erzz=era RIX=RIc RIY=RIb RIZ=RIa erix=eric eriy=erib eriz=eria xx1=c1 yy1=b1 zz1=a1 erxx1=erc1 eryy1=erb1 erzz1=era1 RIX1=RIc1 RIY1=RIb1 RIZ1=RIa1 erix1=eric1 eriy1=erib1 eriz1=eria1 iflagz=1 iflagy=2 else xx=a yy=b zz=c erxx=era eryy=erb erzz=erc RIX=RIa RIY=RIb RIZ=RIc erix=eria eriy=erib eriz=eric xx1=a1 yy1=b1 zz1=c1 erxx1=era1 eryy1=erb1 erzz1=erc1 RIX1=RIa1 RIY1=RIb1 RIZ1=RIc1 erix1=eria1 eriy1=erib1 eriz1=eric1 iflagz=3 iflagy=2 endif c CALL SUBSTS(ZS,YS,erzs,erys,z2,y2) <--- SUBSTS c if(zs.ne.0.d0)then ERCOSZ=DSQRT((0.0015D0/ZS)**2+ERZS**2) else ercosz=0.d0 endif if(ys.ne.0.d0)then ERCOSY=DSQRT((0.0015D0/YS)**2+ERYS**2) else ercosy=0.d0 endif iyy=' ' izz=' ' if(y2.lt.0.d0)iyy='*i' if(z2.lt.0.d0)izz='*i' c WRITE(3,10) WRITE(3,33)' X, Y, Z = ',XX,YY,ZZ write(3,33)' dX, dY, dZ = ',erxx,eryy,erzz WRITE(3,33)' IX, IY, IZ = ',RIx,RIy,RIz WRITE(3,33)'dIX,dIY,dIZ = ',erIx,eRIy,eRIz WRITE(3,35)'Mass = ',TM WRITE(3,15) WRITE(3,33)' X, Y, Z = ',xx1,yy1,zz1 write(3,33)' dX, dY, dZ = ',erxx1,eryy1,erzz1 WRITE(3,33)' IX, IY, IZ = ',RIx1,RIy1,RIz1 WRITE(3,33)'dIX,dIY,dIZ = ',eRIx1,eRIy1,eRIz1 WRITE(3,35)'Mass change = ',DM WRITE(3,35)'Total mass = ',TM+DM WRITE(3,43) c WRITE(*,334)' SYMMETRIC: ',AX(iabs(IFLAgz)),zs,izz ,erzs, * ax(iabs(iflagy)),Ys,iyy ,erys WRITE(*,334)'+Costain err.',AX(iabs(IFLAgz)),zs,izz ,ercosz, * ax(iabs(iflagy)),Ys,iyy ,ercosy WRITE(3,334)' SYMMETRIC: ',AX(iabs(IFLAgz)),zs,izz ,erzs, * ax(iabs(iflagy)),Ys,iyy ,erys WRITE(3,334)'+Costain err.',AX(iabs(IFLAgz)),zs,izz ,ercosz, * ax(iabs(iflagy)),Ys,iyy ,ercosy write(3,'(1x/1x,78(1h-))') goto 124 endif C C...Planar calculation C CALL SUBSTP(XP,YP,erxp,eryp,x2,y2) <--- SUBSTP ERCOSX=DSQRT((0.0015D0/XP)**2+ERXP**2) ERCOSY=DSQRT((0.0015D0/YP)**2+ERYP**2) ixxp=' ' iyyp=' ' if(x2.lt.0.d0)ixxp='*i' if(y2.lt.0.d0)iyyp='*i' c if(iflagx.lt.0.or.iflagy.lt.0)then if(iflagx.lt.0)irecal='I'//ax(iabs(iflagx)) if(iflagy.lt.0)irecal='I'//ax(iabs(iflagy)) write(*,336)rix,riy,irecal,rix1,riy1,irecal srix=rix sriy=riy srix1=rix1 sriy1=riy1 endif WRITE(*,334)' PLANAR: ',AX(iabs(IFLAgx)),XP,ixxp,erxp, * ax(iabs(iflagy)),YP,iyyp,eryp WRITE(*,334)'+Costain err.',AX(iabs(IFLAgx)),XP,ixxp,ercosx, * ax(iabs(iflagy)),YP,iyyp,ercosy 334 FORMAT(1X,A,3X,A1,'=',F10.5,a2,'+-',F7.5, * 4X,A1,'=',F10.5,a2,'+-',F7.5) 336 FORMAT(' Modified values used for PLANAR calculation:'/ * ' parent IX,IY =',2f20.8,' <-- ',a,' recalculated'/ * ' isotopicIX,IY =',2f20.8,' <-- ',a,' recalculated'/) C c...General calculation c xx=a yy=b zz=c erxx=era eryy=erb erzz=erc RIX=RIA RIY=RIB RIZ=RIC erix=eria eriy=erib eriz=eric xx1=a1 yy1=b1 zz1=c1 erxx1=era1 eryy1=erb1 erzz1=erc1 RIX1=RIA1 RIY1=RIB1 RIZ1=RIC1 erix1=eria1 eriy1=erib1 eriz1=eric1 type=' ' CALL SUBST(X,Y,Z,erx,ery,erz,x2,y2,z2,R) <--- SUBST ERCX=DSQRT((0.0015D0/X)**2+ERX**2) ERCY=DSQRT((0.0015D0/Y)**2+ERY**2) ERCZ=DSQRT((0.0015D0/Z)**2+ERZ**2) ERR=dsqrt( (X*ercx)**2+(Y*ercy)**2+(Z*ercz)**2 )/R c ixx=' ' iyy=' ' izz=' ' if(x2.lt.0.d0)ixx='*i' if(y2.lt.0.d0)iyy='*i' if(z2.lt.0.d0)izz='*i' WRITE(*,34)' NONPLANAR: ',X,ixx,erx,Y,iyy,ery,type,Z,izz,erz,R 34 FORMAT(/1X,A,' a=',F10.5,a2,'+-',F7.5, * ' b=',F10.5,a2,'+-',F7.5/1x,a, * ' c=',F10.5,a2,'+-',F7.5,28x,'R=',f8.5) WRITE(*,134)'+Costain err.',X,ixx,ercx,Y,iyy,ercy,err, * type,Z,izz,ercz 134 FORMAT(1X,A,' a=',F10.5,a2,'+-',F7.5, * ' b=',F10.5,a2,'+-',F7.5,' +-',F7.5/1x,a, * ' c=',F10.5,a2,'+-',F7.5) c WRITE(3,10) WRITE(3,33)' X, Y, Z = ',XX,YY,ZZ write(3,33)' dX, dY, dZ = ',erxx,eryy,erzz WRITE(3,33)' IX, IY, IZ = ',RIx,RIy,RIz WRITE(3,33)'dIX,dIY,dIZ = ',erIx,eRIy,eRIz c WRITE(3,33)'PA,PB,PC = ',PA,PB,PC c WRITE(3,33)'In.Defect= ',defect c WRITE(3,33)'Kappa = ',rkappa c WRITE(3,33)'A+B /C = ',APLUSB WRITE(3,35)'Mass = ',TM WRITE(3,15) WRITE(3,33)' X, Y, Z = ',xx1,yy1,zz1 write(3,33)' dX, dY, dZ = ',erxx1,eryy1,erzz1 WRITE(3,33)' IX, IY, IZ = ',RIx1,RIy1,RIz1 WRITE(3,33)'dIX,dIY,dIZ = ',eRIx1,eRIy1,eRIz1 c WRITE(3,33)'PA,PB,PC = ',PA1,PB1,PC1 c WRITE(3,33)'In.Defect= ',defec1 c WRITE(3,33)'Kappa = ',rkapp1 c WRITE(3,33)'A+B /C = ',APLUS1 WRITE(3,35)'Mass change = ',DM WRITE(3,35)'Total mass = ',TM+DM WRITE(3,43) c if(iflagx.lt.0.or.iflagy.lt.0) * write(3,336)srix,sriy,irecal,srix1,sriy1,irecal WRITE(3,334)' PLANAR: ',AX(iabs(IFLAgx)),XP,ixxp,erxp, * ax(iabs(iflagy)),YP,iyyp,eryp WRITE(3,334)'+Costain err.',AX(iabs(IFLAgx)),XP,ixxp,ercosx, * ax(iabs(iflagy)),YP,iyyp,ercosy WRITE(3,34) ' NONPLANAR: ',X,ixx,erx,Y,iyy,ery,type, * Z,izz,erz,R WRITE(3,134)'+Costain err.',X,ixx,ercx,Y,iyy,ercy,err, * type,Z,izz,ercz c write(3,'(1x/1x,78(1h-))') C GOTO 124 c 99 close(3) write(*,'(1x/'' Output has been written to KRA.OUT''/)') STOP c 2000 write(*,2001)errmes 2001 Format(1x//'**** INPUT ERROR: ',a//) c stop END c c------------------------------------------------------------------------ c SUBROUTINE SUBST(X,Y,Z,erx,ery,erz,x2,y2,z2,R) c C...ASYMMETRIC TOP - general case c [Gordy&Cook, eq. 13.71-13.74 p.662] c IMPLICIT REAL*8(A-Z) COMMON IX,IY,IZ,IX1,IY1,IZ1,RM, * erix,eriy,eriz,erix1,eriy1,eriz1 c DIX=IX1-IX DIY=IY1-IY DIZ=IZ1-IZ XY=IX-IY XZ=IX-IZ YZ=IY-IZ DPX=0.5D0*(DIY+DIZ-DIX) DPY=0.5D0*(DIX+DIZ-DIY) DPZ=0.5D0*(DIX+DIY-DIZ) c WRITE(*,33)'DIX,DIY,DIZ = ',DIX,DIY,DIZ C WRITE(*,33)'XY,XZ,YZ = ',XY,XZ,YZ c WRITE(*,33)'DPX,DPY,DPZ = ',DPX,DPY,DPZ 33 FORMAT(1X,A,3F20.5) c X2=DPX*(1+DPY/XY)*(1+DPZ/XZ)/RM X=DSQRT(DABS(X2)) c y2=DPY*(1+DPZ/YZ)*(1-DPX/XY)/RM Y=DSQRT(DABS(y2)) c z2=DPZ*(1-DPX/XZ)*(1-DPY/YZ)/RM Z=DSQRT(DABS(z2)) c A=(DIX+DIY+DIZ)/(2.D0*RM) IF(A.LT.0.D0)WRITE(*,33)' WARNING (negative value) ' R=DSQRT(DABS(A)) c c c...Errors in Kraitchman quantities Delta(P.x), Delta(P.y) and Delta(P.z) c erdpx=0.5d0*dsqrt(erix**2+eriy**2+eriz**2+erix1**2+eriy1**2+ $ eriz1**2) erdpy=erdpx erdpz=erdpx c c 2 2 2 c....Errors in squares of Kraitchman coordinates ie. in X , Y and Z c c The following expressions for the error in the square of the Kraitchman c coordinate have been evaluated with DERIVE by starting from the two c expressions commented out for each error c c X2:=dpx*(dpz+ix-iz)*(dpy+ix-iy)/(rm*(ix-iy)*(ix-iz)) c c DIF(X2,dpx)**2*erdpx**2+DIF(X2,dpy)**2*erdpy**2+DIF(X2,dpz)**2*er c $dpz**2+DIF(X2,ix)**2*erix**2+DIF(X2,iy)**2*eriy**2+DIF(X2,iz)**2*e c $riz**2 c erx2= $dpx**2*(dpy**2*(dpz**2*(erix**2*(2*ix-iy-iz)**2+eriy**2*(ix-iz)**2 $+eriz**2*(ix-iy)**2)+2*dpz*(ix-iz)**2*(erix**2*(2*ix-iy-iz)+eriy** $2*(ix-iz))+(ix-iz)**2*(erdpz**2*(ix-iy)**2+(ix-iz)**2*(erix**2+eri $y**2)))+2*dpy*(ix-iy)**2*(dpz**2*(erix**2*(2*ix-iy-iz)+eriz**2*(ix $-iy))+dpz*erix**2*(ix-iz)**2+erdpz**2*(ix-iz)**2*(ix-iy))+(ix-iy)* $*2*(dpz**2*(erdpy**2*(ix-iz)**2+(ix-iy)**2*(erix**2+eriz**2))+2*dp $z*erdpy**2*(ix-iz)**3+(ix-iz)**2*(erdpy**2*(ix**2-2*ix*iz+iz**2)+e $rdpz**2*(ix**2-2*ix*iy+iy**2))))/(rm**2*(ix-iy)**4*(ix-iz)**4)+erd $px**2*(dpz+ix-iz)**2*(dpy+ix-iy)**2/(rm**2*(ix-iy)**2*(ix-iz)**2) erx2=dsqrt(erx2) c c Y2:=dpy/rm*(1+dpz/(iy-iz))*(1+dpx/(iy-ix)) c c DIF(Y2,dpx)**2*erdpx**2+DIF(Y2,dpy)**2*erdpy**2+DIF(Y2,dpz)**2*er c $dpz**2+DIF(Y2,ix)**2*erix**2+DIF(Y2,iy)**2*eriy**2+DIF(Y2,iz)**2*e c $riz**2 c ery2= $(dpx**2*(dpy**2*(dpz**2*(erix**2*(iy-iz)**2+eriy**2*(ix-2*iy+iz)** $2+eriz**2*(ix-iy)**2)+2*dpz*(iy-iz)**2*(erix**2*(iy-iz)+eriy**2*(- $ix+2*iy-iz))+(iy-iz)**2*(erdpz**2*(ix-iy)**2+(iy-iz)**2*(erix**2+e $riy**2)))+erdpy**2*(iy-iz)**2*(ix-iy)**2*(dpz+iy-iz)**2)-2*dpx*(ix $-iy)**2*(dpy**2*(dpz**2*(eriy**2*(ix-2*iy+iz)+eriz**2*(ix-iy))-dpz $*eriy**2*(iy-iz)**2+erdpz**2*(iy-iz)**2*(ix-iy))+erdpy**2*(iy-iz)* $*2*(ix-iy)*(dpz+iy-iz)**2)+(ix-iy)**2*(dpy**2*(dpz**2*(erdpx**2*( $iy-iz)**2+(ix-iy)**2*(eriy**2+eriz**2))+2*dpz*erdpx**2*(iy-iz)**3 $+(iy-iz)**2*(erdpz**2*(ix-iy)**2+erdpx**2*(iy**2-2*iy*iz+iz**2))) $+erdpy**2*(iy-iz)**2*(ix-iy)**2*(dpz+iy-iz)**2))/(rm**2*(ix-iy)**4 $*(iy-iz)**4) ery2=dsqrt(ery2) c c Z2:=dpz/rm*(1+dpx/(iz-ix))*(1+dpy/(iz-iy)) c c DIF(Z2,dpx)**2*erdpx**2+DIF(Z2,dpy)**2*erdpy**2+DIF(Z2,dpz)**2*er c $dpz**2+DIF(Z2,ix)**2*erix**2+DIF(Z2,iy)**2*eriy**2+DIF(Z2,iz)**2*e c $riz**2 c erz2= $(dpx**2*(dpy**2*(dpz**2*(erix**2*(iy-iz)**2+eriy**2*(ix-iz)**2+eri $z**2*(ix+iy-2*iz)**2)+erdpz**2*(iy-iz)**2*(ix-iz)**2)-2*dpy*(iy-iz $)**2*(dpz**2*(erix**2*(iy-iz)+eriz**2*(ix+iy-2*iz))+erdpz**2*(iy-i $z)*(ix-iz)**2)+(iy-iz)**2*(dpz**2*(erdpy**2*(ix-iz)**2+(iy-iz)**2* $(erix**2+eriz**2))+erdpz**2*(iy-iz)**2*(ix-iz)**2))-2*dpx*(ix-iz)* $*2*(dpy**2*(dpz**2*(eriy**2*(ix-iz)+eriz**2*(ix+iy-2*iz))+erdpz**2 $*(iy-iz)**2*(ix-iz))-dpy*(iy-iz)**2*(dpz**2*eriz**2+2*erdpz**2*(iy $-iz)*(ix-iz))+(iy-iz)**2*(ix-iz)*(dpz**2*erdpy**2+erdpz**2*(iy-iz) $**2))+(ix-iz)**2*(dpy**2*(dpz**2*(erdpx**2*(iy-iz)**2+(ix-iz)**2* $(eriy**2+eriz**2))+erdpz**2*(iy-iz)**2*(ix-iz)**2)-2*dpy*(iy-iz)** $3*(dpz**2*erdpx**2+erdpz**2*(ix-iz)**2)+(iy-iz)**2*(dpz**2*(erdpy $**2*(ix-iz)**2+erdpx**2*(iy-iz)**2)+erdpz**2*(iy-iz)**2*(ix-iz)** $2)))/(rm**2*(ix-iz)**4*(iy-iz)**4) erz2=dsqrt(erz2) c c...Errors in Kraitchman coordinates X, Y and Z c erx=0.5d0*erx2/x ery=0.5d0*ery2/y erz=0.5d0*erz2/z c c RETURN END c c------------------------------------------------------------------------ c SUBROUTINE SUBSTS(Z,Y,erz,ery,z2,y2) c C...SYMMETRIC TOP [on axis: Gordy&Cook, eq. 13.46, p.658] c [off-axis: Gordy&Cook, eq. 13.56,13.57, p.660] c IMPLICIT REAL*8(A-Z) COMMON IX,IY,IZ,IX1,IY1,IZ1,RM, * erix,eriy,eriz,erix1,eriy1,eriz1 c c...substitution on symmetry axis c if((ix1.eq.iy1).or.(iy1.eq.iz1))then DIy=Iy1-Iy z2=diy/rm z=dsqrt(abs(z2)) ERZ2=DSQRT(ERIY1**2+ERIY**2) erz=erz2/(2.d0*z) y=0 y2=0 ery=0 return endif c c...substitution off the symmetry axis c y2=((ix1-iy1)*(ix1-iz1)) / (rm * (ix1-iy1-iz1+iy)) z2=((iy1-iy )*(iz1-iy )) / (rm * (iy1+iz1-ix1-iy)) z=dsqrt(abs(z2)) y=dsqrt(abs(y2)) C ERY2=DSQRT( $erix1**2*(ix1**2+2*ix1*(iy-iy1-iz1)+iy*(-iy1-iz1)+iy1**2+iy1*iz1+i $z1**2)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy**2*(ix1-iz1)**2*(ix1-iy $1)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy1**2*(iy-iz1)**2*(ix1-iz1)** $2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriz1**2*(iy-iy1)**2*(ix1-iy1)**2/(r $m**2*(ix1+iy-iy1-iz1)**4) ) ery=erz2/(2.d0*y) ERZ2=DSQRT( $erix1**2*(iy-iz1)**2*(iy-iy1)**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy* $*2*(ix1*(2*iy-iy1-iz1)+iy**2-2*iy*(iy1+iz1)+iy1**2+iy1*iz1+iz1**2) $**2/(rm**2*(ix1+iy-iy1-iz1)**4)+eriy1**2*(iy-iz1)**2*(ix1-iz1)**2/ $(rm**2*(ix1+iy-iy1-iz1)**4)+eriz1**2*(iy-iy1)**2*(ix1-iy1)**2/(rm* $*2*(ix1+iy-iy1-iz1)**4) ) erz=erz2/(2.d0*z) c c RETURN END c c------------------------------------------------------------------------ c SUBROUTINE SUBSTP(X,Y,erx,ery,x2,y2) c C..ASYMMETRIC TOP P L A N A R C [Gordy&Cook, eq. 13.64,13.65 p.661] C C All input parameters are passed through the COMMON blocks c c iflagx, iflagy - inertial coordinates to be treated as x and y c 1,2,3 for a,b,c. A negative value can be specified for c one of iflagx,iflagy, in which case the necessary quantities c for this axis are evaluated from moments for the remaining c two axes using the planarity relations: c I.b+ I.a= I.c c DI.b+DI.a=DI.c c as defined for the case in question by ifirst,isecnd, ie. if c iflagx is negative then c I.iflagx = I.ifirst + I.isecnd c c X, Y - r.s coordinates from planar equations c ERX, ERY - errors on r.s coordinates from error propagation c X2, Y2 - squares of r.s coordinates as calculated with the equations, c these can be negative which can be used for marking imaginary c coordinates c IMPLICIT REAL*8(A-Z) integer iflagx,iflagy,ifirst,isecnd,i1,i2 COMMON IX,IY,IZ,IX1,IY1,IZ1,RM , * erix,eriy,eriz,erix1,eriy1,eriz1 common /moment/parent(3,2),isot(3,2),iflagx,iflagy,ifirst,isecnd c if(iflagx.gt.0)then ix=parent(iflagx,1) ix1=isot(iflagx,1) DIX=IX1-IX ERIX=PARENT(IFLAGX,2) ERDIX=DSQRT(PARENT(IFLAGX,2)**2+ISOT(IFLAGX,2)**2) else i1=iabs(ifirst) i2=iabs(isecnd) s1=isign(1,ifirst) s2=isign(1,isecnd) dix=s1*(isot(i1,1)-parent(i1,1)) * +s2*(isot(i2,1)-parent(i2,1)) ix=s1*parent(i1,1)+s2*parent(i2,1) ERDIX=DSQRT(ISOT(i1,2)**2+PARENT(i1,2)**2 * +ISOT(i2,2)**2+PARENT(i2,2)**2) ERIX=DSQRT(PARENT(i1,2)**2+PARENT(i2,2)**2) endif c if(iflagy.gt.0)then iy=parent(iflagy,1) iy1=isot(iflagy,1) DIY=IY1-IY ERIY=PARENT(IFLAGY,2) ERDIY=DSQRT(PARENT(IFLAGY,2)**2+ISOT(IFLAGY,2)**2) else i1=iabs(ifirst) i2=iabs(isecnd) s1=isign(1,ifirst) s2=isign(1,isecnd) diy=s1*(isot(i1,1)-parent(i1,1)) * +s2*(isot(i2,1)-parent(i2,1)) iy=s1*parent(i1,1)+s2*parent(i2,1) ERDIy=DSQRT(ISOT(i1,2)**2+PARENT(i1,2)**2 * +ISOT(i2,2)**2+PARENT(i2,2)**2) ERIy=DSQRT(PARENT(i1,2)**2+PARENT(i2,2)**2) endif XY=IX-IY ERXY=DSQRT(ERIX**2+ERIY**2) c WRITE(*,33)'DIX,DIY,XY = ',DIX,DIY,XY c WRITE(3,33)'DIX,DIY,XY = ',DIX,DIY,XY c33 FORMAT(1X,A,3F20.5) C x2=DIY*(1.d0+DIX/XY)/RM X=DSQRT(DABS(x2)) y2=DIX*(1.d0-DIY/XY)/RM Y=DSQRT(DABS(y2)) c erx2= *DSQRT(dix**2*(diy**2*erxy**2+erdiy**2*xy**2)+2*dix*erdiy**2*xy**3 * +diy**2*erdix**2*xy**2+erdiy**2*xy**4)/(xy**2*ABS(rm)) erx=erx2/(2.d0*x) ery2= *DSQRT(dix**2*(diy**2*erxy**2+erdiy**2*xy**2)+diy**2*erdix**2*xy**2 * -2*diy*erdix**2*xy**3+erdix**2*xy**4)/(xy**2*ABS(rm)) ery=ery2/(2.d0*y) c RETURN END c c------------------------------------------------------------------------ c subroutine test(iodev,ierr) character line*80 c c...test whether the next line in currently open file on IOUNIT is a comment c line (IERR=1), otherwise IERR=0 c read(iodev,'(a)',err=3,end=3)line if(line(1:1).eq.'c'.or.line(1:1).eq.'C'.or.line(1:1).eq.'!'. *or.line(1:1).eq.'#'.or.line(1:1).eq.'%')then ierr=1 write(3,'(1x,a)')line(2:) else backspace(iodev) ierr=0 endif return c 3 ierr=2 return end c c------------------------------------------------------------------------ c------------------------------------------------------------------------