C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C CONVC - PROGRAM FOR INTERCONVERSION BETWEEN VARIOUS TYPES OF C ASYMMETRIC TOP CONSTANTS C C WATSON S-reduced --> WATSON A-reduced (sextic level + errors) C WATSON A-reduced --> WATSON S-reduced (sextic level + errors) C vanEIJCK/TYPKE S --> WATSON A-reduced (quartics + errors) C KIRCHOFF --> WATSON (quartics only) C WATSON A, I.r --> WATSON A, III.r (quartics only) C WATSON A, I.r --> WATSON A, II.r (quartics only) C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C - Constants are to be input in the .CON format used by ASROT containing c the number of constants in the first line, followed by lines with c constant name, value and error c c - output is produced to: c a/ screen c b/ file CONSTC.OUT which echoes the screen c c/ file CONSTC.CON which is in the .CON format c (both disk files are written to in the APPEND mode) c c - note that conversion will tend to give greater errors than a refit c of the data C c for formulae see: c c 1/ Watson, "Vibrational Spectra and structure", Ed. J.Durig, Vol.6 (1977) c 2/ Yamada, Winnewisser - Z.Naturforsch. 31a,131(1975) c 3/ (Notebook 2,p.6-10) c C ver. 19.IV.2022 ----- Zbigniew KISIEL ----- C __________________________________________________ C | Institute of Physics, Polish Academy of Sciences | C | Al.Lotnikow 32/46, Warszawa, POLAND | C | kisiel@ifpan.edu.pl | C_________________________/-------------------------------------------------- c C xx.10.94: Version with five conversion options and A->S errors c 27.09.98: Elimination of SIGMA inaccuracy in A->S conversion c Addition of S->A errors c 28.09.98: Debugging of h2 calculation c VanEijck/Typke errors c 10.12.98: Kirchhoff->Watson A,Ir c 11.12.98: output mods. c 19.04.22: overhaul and updates c C------------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z) c IMPLICIT INTEGER*2 (I-N) COMMON /DAT/AA(15),AS(15),AERROR(15),SERROR(15) REAL*8 AINP(15),errinp(15),AK(8),ERRKIR(8),VE(15),verror(15), * AAII(15),AAIII(15) integer units(15) CHARACTER FILNAM*30,labin*15,labout*15 CHARACTER*6 WATSA(15),WATSS(15),KIRCH(8),VANEIC(8) CHARACTER*7 UNIT,NAMUN(6) CHARACTER*40 OUTCON(2) CHARACTER*27 CONVAL,ERVAL,CONOUT,EROUT character lintest*200 C EQUIVALENCE (AA(1),A),(AA(2),B),(AA(3),C),(AA(4),DJ), * (AA(5),DJK),(AA(6),DK),(AA(7),DDJ),(AA(8),DDK), * (AA(9),HJ),(AA(10),HJK),(AA(11),HKJ),(AA(12),HK), * (AA(13),HHJ),(AA(14),HHJK),(AA(15),HHK) EQUIVALENCE (AS(1),SA),(AS(2),SB),(AS(3),SC),(AS(4),SDJ), * (AS(5),SDJK),(AS(6),SDK),(AS(7),D1),(AS(8),D2), * (AS(9),SHJ),(AS(10),SHJK),(AS(11),SHKJ),(AS(12),SHK), * (AS(13),H1),(AS(14),H2),(AS(15),H3) EQUIVALENCE (AK(1),AKIR),(AK(2),BKIR),(AK(3),CKIR), * (AK(4),TAAAA),(AK(5),TBBBB),(AK(6),TCCCC),(AK(7),T1), * (AK(8),T2) EQUIVALENCE (errKir(1),eAKIR),(errKir(2),eBKIR), * (errKir(3),eCKIR), * (errKir(4),eTAAAA),(errKir(5),eTBBBB),(errKir(6),eTCCCC), * (errKir(7),eT1),(errKir(8),eT2) EQUIVALENCE (ve(1),vea),(ve(2),veB),(ve(3),veC),(ve(4),vedJ), * (ve(5),vedJK),(ve(6),veDK),(ve(7),veddj),(ve(8),VER6) EQUIVALENCE (verror(1),evea),(verror(2),eveB),(verror(3),eveC), * (verror(4),evedJ),(verror(5),evedJK),(verror(6),eveDK), * (verror(7),eveddj),(verror(8),eVER6) EQUIVALENCE (AERROR(1),EA),(AERROR(2),EB),(AERROR(3),EC), * (AERROR(4),EDJ),(AERROR(5),EDJK),(AERROR(6),EDK), * (AERROR(7),EDDJ),(AERROR(8),EDDK),(AERROR(9),EHJ), * (AERROR(10),EHJK),(AERROR(11),EHKJ),(AERROR(12),EHK), * (AERROR(13),EHHJ),(AERROR(14),EHHJK),(AERROR(15),EHHK) DATA WATSA/'A =','B =','C =',' DJ =',' DJK =', 1 ' DK =',' dJ =',' dK =','HJ =','HJK =','HKJ =', 1 'HK =','hj =','hjk =','hk ='/ DATA WATSS/'A =','B =','C =',' DJ =',' DJK =', 1 ' DK =',' d1 =',' d2 =','HJ =','HJK =','HKJ =', 1 'HK =','h1 =','h2 =','h3 ='/ DATA KIRCH/'A =','B =','C =','Taaaa=','Tbbbb=', 1 'Tcccc=','Tau1 =','Tau2 ='/ DATA VANEIC/'A'' =','B'' =','C'' =',' DJ'' =',' DJK''=', 1 ' DK'' =',' dJ'' =',' R6'' '/ C C...UNITS defines the units in which various constants are to be printed C out, NAMUN contains the names of these units DATA NAMUN/'MHz ','kHz ',' Hz ','mHz ','microHz', * 'nHz '/ DATA UNITS/1,1,1, 2,2,2,2,2, 3,3,3,3,3,3,3/ c CON=505379.0089D0 ! CODATA 2014 nerd=2 C WRITE(*,1) 1 FORMAT(1X//' ',76(1H_)/' |',T79,'|'/ * ' | CONVC - INTERCONVERSION BETWEEN VARIOUS TYPES OF ASYM', * 'METRIC',T79,'|'/ * ' | ROTOR CONSTANTS', * T79,'|'/' |',76(1H_),'|'/ * ' version 19.IV.2022',T64,'Zbigniew KISIEL'//) c 2908 WRITE(*,6855) 6855 FORMAT(/' NAME OF FILE CONTAINING CONSTANTS: '\) 2909 READ(*,'(A)',ERR=2908)FILNAM OPEN(3,FILE=FILNAM,ERR=2908,status='old') READ(3,6866,ERR=2909)NCONST 6866 FORMAT(7X,I5) if(nconst.lt.1.or.nconst.gt.36)goto 2908 IF(NCONST.GT.15)NCONST=15 DO 6856 I=1,NCONST 2911 READ(3,6867,ERR=2910)AINP(I),errinp(i) 6867 FORMAT(7x,2F30.19) GOTO 6856 2910 WRITE(*,'(1X/1X,A,''ERROR IN CONSTANT '',I1)')CHAR(7),I GOTO 2911 6856 CONTINUE CLOSE(3) OPEN(4,FILE='CONVC.OUT',ACCESS='APPEND') OPEN(2,FILE='CONVC.CON',ACCESS='APPEND') C IF(NCONST.GE.15)GOTO 11 DO 3 I=NCONST+1,15 AINP(I)=0.D0 errinp(i)=0.d0 3 CONTINUE C C 11 WRITE(*,10) 10 FORMAT(1X//' SELECT CONVERSION OPTION:'/ * 25X, '0 = E X I T'/ * 25X, '1 = Watson S --> Watson A (I.r)'/ * 25X, '2 = Watson A --> Watson S (I.r)'/ * 25X, '3 = vanEIJCK,TYPKE --> Watson A'/ * 25X, '4 = Kirchhoff --> Watson A'/ * 25X, '5 = Watson A, I.r --> Watson A, III.r'/ * 25X, '6 = Watson A, I.r --> Watson A, II.r'/ * 60X,'..... ',$) READ(*,'(I5)',ERR=11)ICONV IF(ICONV.EQ.0)THEN CLOSE(4) STOP ENDIF IF(ICONV.LT.0.OR.ICONV.GT.6)GOTO 11 GOTO(51,151,52,53,54,55),ICONV C c- - - - - - - - - - - - - - - - - - - - - - - - - - C C...Watson S --> Watson A conversion (with formulae from Table 8.16, p.335 C of 2nd Ed. of Gordy and Cook - reprinted from JMS 65,123(1977)) I.r C 51 labin='Watson S,I.r' labout='Watson A,I.r' DO 12 I=1,15 AS(I)=AINP(I) serror(i)=errinp(i) 12 CONTINUE C SIGMA=(2.D0*SA-(SB+SC) )/(SB-SC) A=SA+10.D0*D2 B=SB-4.D0*(2.d0*SIGMA+1.D0)*D2 C=SC+4.D0*(2.D0*SIGMA-1.D0)*D2 c DJ=SDJ-2.D0*D2 DJK=SDJK +12.D0*D2 DK=SDK -10.D0*D2 DDJ=-D1 DDK=-4.D0*D2*SIGMA c HJ= SHJ +2.D0*H2 HJK=SHJK -12.D0*H2+16.D0*SIGMA*H3 HKJ=SHKJ +10.D0*H2-(160.D0*SIGMA/3.D0)*H3 HK= SHK +(112.0D0*SIGMA/3.D0)*H3 HHJ=H1+H3 HHJK=4.0D0*SIGMA*H2-10.D0*H3-8.D0*D2*(SDJK +2.D0*SIGMA*D1 * +4.D0*D2)/(SB-SC) HHK=(32.D0*SIGMA**2/3.D0+9.D0)*H3-16.D0*D2*(SDK * +2.D0*(SIGMA**2-2.D0)*D2)/(SB-SC) C call errsta <----- c WRITE(*,14) 14 FORMAT(1X/) nout=nconst write(2,18)nout 18 format(' NCON =',i5) write(4,190)labin,labout 190 format(1x/10x,a,t33,'--->',t40,a//) DO 15 I=1,nout WRITE(*,16)WATSS(I),AS(I),WATSA(I),AA(I) WRITE(4,16)WATSS(I),AS(I),WATSA(I),AA(I) write(2,17)watsa(i),aa(i),aerror(i) 15 CONTINUE write(2,19) labout,labin 19 format(' |'/' -- ',A,' calculated from ',A//) 16 FORMAT(1X,A6,F25.15,8X,A6,F25.15) 17 format(1x,a6,2F30.19) c WRITE(4,105) 105 FORMAT(1X/1X,78(1H-)/) GOTO 11 C c- - - - - - - - - - - - - - - - - - - - - - - - - - C C...Watson A --> Watson S (I.r), inversion of the formulae used for S->A C conversion C 151 labin='Watson A,I.r' labout='Watson S,I.r' DO 152 I=1,15 AA(I)=AINP(I) aerror(i)=errinp(i) 152 CONTINUE C c...SIGMA is defined in terms of (initially unknown) S-reduction rotational c constants, which involve a d2 correction that is directly c dependent on SIGMA. c To obtain full A<->S reversability two initial conversion iterations c are made in order to obtain reliable SIGMA c SIGMA=(2.D0*A-(B+C) )/(B-C) D2=-DDK /(4.D0*SIGMA) SA=A-10.D0*D2 SB=B+4.D0*(2.D0*SIGMA+1.D0)*D2 SC=C-4.D0*(2.D0*SIGMA-1.D0)*D2 SIGMA=(2.D0*SA-(SB+SC) )/(SB-SC) D2=-DDK /(4.D0*SIGMA) SA=A-10.D0*D2 SB=B+4.D0*(2.D0*SIGMA+1.D0)*D2 SC=C-4.D0*(2.D0*SIGMA-1.D0)*D2 SIGMA=(2.D0*SA-(SB+SC) )/(SB-SC) C D1=-DDJ D2=-DDK /(4.D0*SIGMA) SDJ= DJ +2.D0*D2 SDJK=DJK-12.D0*D2 SDK= DK +10.D0*D2 SA=A-10.D0*D2 SB=B+4.D0*(2.D0*SIGMA+1.D0)*D2 SC=C-4.D0*(2.D0*SIGMA-1.D0)*D2 C H3=(hhk+16.d0*d2*(sdk+2.d0*(sigma**2-2.d0)*d2)/(sb-sc))/ * (32.d0*sigma**2/3.d0+9.d0) H2=(hhjk+10.d0*h3+8.d0*d2*(sdjk+2.d0*sigma*d1+4.d0*d2)/ * (sb-sc))/(4.d0*sigma) H1=hhj-h3 SHJ=hj-2.d0*h2 SHJK=hjk+12.d0*h2-16.d0*sigma*h3 SHKJ=hkj-10.d0*h2+(160.d0*sigma/3.d0)*h3 SHK=hk-(112.d0*sigma/3.d0)*h3 C CALL ERRATS <----- C WRITE(*,14) nout=nconst write(2,18)nout write(4,190)labin,labout DO 154 N=1,nout WRITE(OUTCON(1),'(38(1H ))') WRITE(OUTCON(2),'(38(1H ))') UNIT=NAMUN(UNITS(N)) ASCLD=AA(N)*10.D0**(3*(UNITS(N)-1)) WRITE(CONVAL,'(F27.16)')ASCLD IF(AERROR(N).NE.0.0D0)THEN ESCLD=AERROR(N)*10.D0**(3*(UNITS(N)-1)) WRITE(ERVAL,'(F27.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR,nerd) <----- c c string length check: 26 = 40 - 4 (name) - 7 (unit) - 3 (/ and () in format 411) c if(ndcon+nderor.gt.26)then ! debug WRITE(lintest,518)CONVAL(1:NDCON),ERVAL(1:NDEROR) ! debug 518 format(a,'(',a,')') ! debug write(*,'(1x,a/1x,a)') ! debug * 'WARNING --> string length overflow'// ! debug * ' on an unusual value of parameter:', ! debug * lintest(1:len_trim(lintest)) ! debug write(*,'(1x,a//1x,a,3x)',advance='NO') ! debug * 'Problems expected - untested condition !', ! debug * 'Press ENTER to continue' ! debug read(*,*,err=9999)idummy ! debug 9999 continue ! debug erval='0' ! debug nderor=1 ! debug endif ! debug c WRITE(OUTCON(1),411)WATSA(N)(1:4),UNIT,CONVAL(1:NDCON), * ERVAL(1:NDEROR) ELSE WRITE(OUTCON(1),1005)WATSA(N)(1:4),UNIT,CONVAL(1:27) ENDIF C ASCLD=AS(N)*10.D0**(3*(UNITS(N)-1)) WRITE(CONVAL,'(F27.16)')ASCLD IF(SERROR(N).NE.0.D0)THEN ESCLD=SERROR(N)*10.D0**(3*(UNITS(N)-1)) WRITE(ERVAL,'(F27.16)')ESCLD CALL CONFOR(CONVAL,ERVAL,NDCON,NDEROR,nerd) <----- c c string length check: 26 = 40 - 4 (name) - 7 (unit) - 3 (/ and () in format 411) c if(ndcon+nderor.gt.26)then ! debug WRITE(lintest,518)CONVAL(1:NDCON),ERVAL(1:NDEROR) ! debug write(*,'(1x,a/1x,a)') ! debug * 'WARNING --> string length overflow'// ! debug * ' on an unusual value of parameter:', ! debug * lintest(1:len_trim(lintest)) ! debug write(*,'(1x,a//1x,a,3x)',advance='NO') ! debug * 'Problems expected - untested condition !', ! debug * 'Press ENTER to continue' ! debug read(*,*,err=9998)idummy ! debug 9998 continue ! debug erval='0' ! debug nderor=1 ! debug endif ! debug c WRITE(OUTCON(2),411)WATSS(N)(1:4),UNIT,CONVAL(1:NDCON), * ERVAL(1:NDEROR) ELSE WRITE(OUTCON(2),1005)WATSS(N)(1:4),UNIT,CONVAL(1:27) ENDIF C WRITE(*,1108)OUTCON(1),OUTCON(2) WRITE(4,1108)OUTCON(1),OUTCON(2) write(2,17)watss(n),as(n),serror(n) 154 CONTINUE c write(2,19)labout,labin 1108 FORMAT(1X,A38,1X,A38) 411 FORMAT(A,'/',A,A,'(',A,') ') 1005 FORMAT(A,'/',A,A) C WRITE(4,105) GOTO 11 c c- - - - - - - - - - - - - - - - - - - - - - - - - - C C...vanEIJCK/TYPKE S-reduction ---> WATSON A-REDUCTION C 52 labin='Van Eijck/Typke' labout='Watson A,I.r' DO 20 I=1,8 VE(I)=AINP(I) verror(i)=errinp(i) 20 CONTINUE C DDKPR=-2.D0*VER6 DDK =-2.D0*DDKPR*(2.D0*VEA-VEB-VEC)/(VEB-VEC) A =VEA-8.D0*DDKPR B =VEB+4.D0*DDKPR+2.D0*DDK C =VEC+4.D0*DDKPR-2.D0*DDK DJ =VEDJ +DDKPR DJK =VEDJK-6.D0*DDKPR DK =VEDK +5.D0*DDKPR DDJ =VEDDJ C EDDKPR=2.D0*EVER6 eddk=dsqrt( $ ddkpr**2*(16.d0*evea**2/(veb-vec)**2+16.d0*eveb**2*(vea-vec)**2 $ /(veb-vec)**4+16.d0*evec**2*(vea-veb)**2/(veb-vec)**4) $ +4.d0*eddkpr**2*(2.d0*vea-veb-vec)**2/(veb-vec)**2 ) ea=dsqrt( 64.d0*eddkpr**2+evea**2 ) eb=dsqrt( 16.d0*eddkpr**2+eveb**2+4.d0*eddk**2 ) ec=dsqrt( 16.d0*eddkpr**2+evec**2+4.d0*eddk**2 ) c edj=dsqrt( eddkpr**2+evedj**2 ) edjk=dsqrt( 36.d0*eddkpr**2+evedjk**2 ) edk=dsqrt( 25.d0*eddkpr**2+evedk**2 ) eddj=eveddj C WRITE(*,14) write(2,18)8 write(4,190)labin,labout DO 21 I=1,8 WRITE(*,16)VANEIC(I),VE(I),WATSA(I),AA(I) WRITE(4,16)VANEIC(I),VE(I),WATSA(I),AA(I) write(2,17)watsa(i),aa(i),aerror(i) 21 CONTINUE write(2,19)labout,labin WRITE(4,105) GOTO 11 c c- - - - - - - - - - - - - - - - - - - - - - - - - - C C....Kirchhoff -> Watson A I.r C C 53 labin='Kirchhoff' labout='Watson A,I.r' IF(NCONST.GT.8)NCONST=8 DO 200 I=1,8 AK(I)=AINP(I) ERRKIR(I)=ERRINP(I) AERROR(I)=0.D0 200 CONTINUE C A =AKIR B =BKIR C =CKIR C DJ= -0.125D0*(TBBBB+TCCCC) DJK= 0.125D0*(3.D0*TBBBB+3.D0*TCCCC-2.D0*T1) DK= -0.25D0 *(TAAAA+TBBBB+TCCCC-T1) DDJ=-0.0625D0*(TBBBB-TCCCC) DDK= 0.125D0/(B-C)*( (B-A)*TBBBB+(C-A)*TCCCC-(B+C)*T1 * +2.D0*(A+B+C)*T2 ) A =AKIR-2.D0*DJ B =BKIR-2.D0*DJ-DJK+2.D0*(DDJ+DDK) C =CKIR-2.D0*DJ-DJK-2.D0*(DDJ+DDK) C edj= 0.125d0*dsqrt(etbbbb**2+etcccc**2) edjk= 0.125d0*dsqrt(9.d0*etbbbb**2+9.d0*etcccc**2+4.d0*et1**2) edk= 0.25d0*dsqrt(etaaaa**2+etbbbb**2+etcccc**2+et1**2) eddj=0.0625d0*dsqrt(etbbbb**2+etcccc**2) eddk= dsqrt( $(a**2*(2*t2-tbbbb-tcccc)**2*(eb**2+ec**2)-2*a*(2*t2-tbbbb-tcccc)*( $2*t1-4*t2-tbbbb-tcccc)*(b*ec**2+c*eb**2)+b**2*(ea**2*(2*t2-tbbbb-t $cccc)**2+ec**2*(2*t1-4*t2-tbbbb-tcccc)**2)-2*b*c*ea**2*(2*t2-tbbbb $-tcccc)**2+c**2*(ea**2*(2*t2-tbbbb-tcccc)**2+eb**2*(2*t1-4*t2-tbbb $b-tcccc)**2))/(64*(b-c)**4)+(a**2*(4*et2**2+etbbbb**2+etcccc**2)+2 $*a*(b*(4*et2**2-etbbbb**2)+c*(4*et2**2-etcccc**2))+b**2*(et1**2+4* $et2**2+etbbbb**2)+2*b*c*(et1**2+4*et2**2)+c**2*(et1**2+4*et2**2+et $cccc**2))/(64*(b-c)**2) ) c ea=dsqrt(eakir**2+4.d0*edj**2) eb=dsqrt(ebkir**2+4.d0*edj**2+edjk**2+4.d0*eddj**2+4.d0*eddk**2) ec=dsqrt(eckir**2+4.d0*edj**2+edjk**2+4.d0*eddj**2+4.d0*eddk**2) c WRITE(*,14) write(2,18)8 write(4,190)labin,labout DO 201 I=1,8 WRITE(*,16)KIRCH(I),AK(I),WATSA(I),AA(I) WRITE(4,16)KIRCH(I),AK(I),WATSA(I),AA(I) write(2,17)watsa(i),aa(i),aerror(i) 201 CONTINUE write(2,19)labout,labin WRITE(4,105) GOTO 11 C goto 101 c c- - - - - - - - - - - - - - - - - - - - - - - - - - C C....Watson A,I.r -> Watson A,III.r c (reverse the signs of resulting delta.J and C delta.K for Watson I.r -> Watson III.l) C This uses formulae from Yamada, Winnewisser - Z.Naturforsch. 31a,131(1975): C C Eq.10: T0 = Q_i D_i where i=I,II,and III (i.e. T0 is axis independent) C Hence T0 = Q_I D_I = Q_III D_III C Leading to D_III = Q_III^(-1) Q_I D_I C where Q_III^(-1) Q_I is given by Eq(17) and Eq.(a8) C and D_I is a column vector [DeltaJ, DeltaJK, DeltaK, deltaJ, deltaK]_I C similarly D_III C 54 labin='Watson A,I.r' labout='Watson A,III.r' DO 102 I=1,15 AA(I)=AINP(I) Aerror(i)=0.d0 102 CONTINUE c aaiii(1)=aa(1) aaiii(2)=aa(2) aaiii(3)=aa(3) aaiii(4)= 1.0d0*dj + 0.5d0*djk + 0.5d0*dk + ddj aaiii(5)= - 0.5d0*djk - 1.5d0*dk - 3.d0*ddj aaiii(6)= dk aaiii(7)= 0.25d0*djk + 0.25d0*dk - 0.5d0*ddj aaiii(8)= -(a-c)/(2.d0*(a-b))*dk + (b-c)/(a-b)*ddk c c...OUT file output c WRITE(*,14) write(4,190)labin,labout DO 100 I=1,8 WRITE(*,16)WATSA(I),AA(I),WATSA(I),AAiii(I) WRITE(4,16)WATSA(I),AA(I),WATSA(I),AAiii(I) 100 CONTINUE 160 FORMAT(1X,A6,F24.14,F19.14,10X,A6,F24.14,F19.14) WRITE(4,105) c c...CON file output c write(2,18)8 DO 1102 I=1,8 write(2,17)watsa(i),aaIII(i),aerror(i) 1102 CONTINUE write(2,19)labout,labin GOTO 11 c c- - - - - - - - - - - - - - - - - - - - - - - - - - C C....Watson A,I.r -> Watson A,II.r C C Also using Yamada, Winnewisser - Z.Naturforsch. 31a,131(1975) (see above): C C D_II = Q_II^(-1) Q_I D_I C where Q_II^(-1) Q_I is given by Eq(a7) C 55 labin='Watson A,I.r' labout='Watson A,II.r' c DO 1105 I=1,15 AA(I)=AINP(I) Aerror(i)=0.d0 1105 CONTINUE c aaii(1)=aa(1) aaii(2)=aa(2) aaii(3)=aa(3) aaii(4)= 1.0d0*dj + 0.5d0*djk + 0.5d0*dk - ddj aaii(5)= - 0.5d0*djk - 1.5d0*dk + 3.d0*ddj aaii(6)= dk aaii(7)= - 0.25d0*djk - 0.25d0*dk - 0.5d0*ddj aaii(8)= (a-b)/(2.d0*(a-c))*dk - (b-c)/(a-c)*ddk c c...OUT file output c WRITE(*,14) write(4,190)labin,labout DO 1100 I=1,8 WRITE(*,16)WATSA(I),AA(I),WATSA(I),AAii(I) WRITE(4,16)WATSA(I),AA(I),WATSA(I),AAii(I) 1100 CONTINUE WRITE(4,105) c c...CON file output c write(2,18)8 DO 1103 I=1,8 write(2,17)watsa(i),aaII(i),aerror(i) 1103 CONTINUE write(2,19)labout,labin GOTO 11 C c- - - - - - - - - - - - - - - - - - - - - - - - - - c 101 STOP END C C------------------------------------------------------------------------- C SUBROUTINE ERRATS C C ERRors on A TO S conversion: C C EVALUATION OF ERRORS IN CALCULATED S-REDUCTION CONSTANTS FROM PROPAGATION C OF ERRORS ON THE INGOING A-REDUCTION CONSTANTS C C formulae from DERIVE and the generating relations for that program have C been left in the comment lines C IMPLICIT REAL*8 (A-H,O-Z) COMMON /DAT/AA(15),AS(15),AERROR(15),SERROR(15) C EQUIVALENCE (AA(1),A),(AA(2),B),(AA(3),C),(AA(4),DJ), * (AA(5),DJK),(AA(6),DK),(AA(7),DDJ),(AA(8),DDK), * (AA(9),HJ),(AA(10),HJK),(AA(11),HKJ),(AA(12),HK), * (AA(13),HHJ),(AA(14),HHJK),(AA(15),HHK) EQUIVALENCE (AS(1),SA),(AS(2),SB),(AS(3),SC),(AS(4),SDJ), * (AS(5),SDJK),(AS(6),SDK),(AS(7),D1),(AS(8),D2), * (AS(9),SHJ),(AS(10),SHJK),(AS(11),SHKJ),(AS(12),SHK), * (AS(13),H1),(AS(14),H2),(AS(15),H3) EQUIVALENCE (AERROR(1),EA),(AERROR(2),EB),(AERROR(3),EC), * (AERROR(4),EDJ),(AERROR(5),EDJK),(AERROR(6),EDK), * (AERROR(7),EDDJ),(AERROR(8),EDDK),(AERROR(9),EHJ), * (AERROR(10),EHJK),(AERROR(11),EHKJ),(AERROR(12),EHK), * (AERROR(13),EHHJ),(AERROR(14),EHHJK),(AERROR(15),EHHK) EQUIVALENCE (SERROR(1),ESA),(SERROR(2),ESB),(SERROR(3),ESC), * (SERROR(4),ESDJ),(SERROR(5),ESDJK),(SERROR(6),ESDK), * (SERROR(7),ED1),(SERROR(8),ED2),(SERROR(9),ESHJ), * (SERROR(10),ESHJK),(SERROR(11),ESHKJ),(SERROR(12),ESHK), * (SERROR(13),EH1),(SERROR(14),EH2),(SERROR(15),EH3) C c c...errors in quadratic and quartic constants c c S:=(2*a-b-c)/(b-c) c D2:=(-ddk/(4*S)) c D1:=(-ddj) c SDJ:=dj+2*D2 c SDJK:=djk-12*D2 c SDK:=dk+10*D2 c SB:=b+4*(2*S+1)*D2 c SC:=c-4*(2*S-1)*D2 c SA:=a-10*D2 c c ea**2*DIF(SA,a)**2+eb**2*DIF(SA,b)**2+ec**2*DIF(SA,c)**2+eddk**2*D c $IF(SA,ddk)**2 c esa= dsqrt( $ea**2*(4*a**2-4*a*(b+c)+b**2+b*(2*c-5*ddk)+c*(c+5*ddk))**2/(2*a-b- $c)**4+25*ddk**2*eb**2*(a-c)**2/(2*a-b-c)**4+25*ddk**2*ec**2*(a-b)* $*2/(2*a-b-c)**4+25*eddk**2*(b-c)**2/(4*(2*a-b-c)**2) ) c c ea**2*DIF(SB,a)**2+eb**2*DIF(SB,b)**2+ec**2*DIF(SB,c)**2+eddk**2*D c $IF(SB,ddk)**2 c esb= dsqrt( $4*ddk**2*ea**2*(b-c)**2/(2*a-b-c)**4+eb**2*(4*a**2-2*a*(2*b+2*c+dd $k)+b**2+2*b*c+c*(c+2*ddk))**2/(2*a-b-c)**4+4*ddk**2*ec**2*(a-b)**2 $/(2*a-b-c)**4+eddk**2*(4*a-b-3*c)**2/(2*a-b-c)**2 ) c c ea**2*DIF(SC,a)**2+eb**2*DIF(SC,b)**2+ec**2*DIF(SC,c)**2+eddk**2*D c $IF(SC,ddk)**2 c esc= dsqrt( $4*ddk**2*ea**2*(b-c)**2/(2*a-b-c)**4+4*ddk**2*eb**2*(a-c)**2/(2*a- $b-c)**4+ec**2*(4*a**2-2*a*(2*b+2*c-ddk)+b**2+2*b*(c-ddk)+c**2)**2/ $(2*a-b-c)**4+eddk**2*(4*a-3*b-c)**2/(2*a-b-c)**2 ) c c ea**2*DIF(SDJ,a)**2+eb**2*DIF(SDJ,b)**2+ec**2*DIF(SDJ,c)**2+eddk** c $2*DIF(SDJ,ddk)**2+edj**2*DIF(SDJ,dj)**2 c esdj= dsqrt( $ddk**2*ea**2*(b-c)**2/(2*a-b-c)**4+ddk**2*eb**2*(a-c)**2/(2*a-b-c) $**4+ddk**2*ec**2*(a-b)**2/(2*a-b-c)**4+eddk**2*(b-c)**2/(4*(2*a-b- $c)**2)+edj**2 ) c c ea**2*DIF(SDJK,a)**2+eb**2*DIF(SDJK,b)**2+ec**2*DIF(SDJK,c)**2+edd c $k**2*DIF(SDJK,ddk)**2+edjk**2*DIF(SDJK,djk)**2 c esdjk= dsqrt( $36*ddk**2*ea**2*(b-c)**2/(2*a-b-c)**4+36*ddk**2*eb**2*(a-c)**2/(2* $a-b-c)**4+36*ddk**2*ec**2*(a-b)**2/(2*a-b-c)**4+9*eddk**2*(b-c)**2 $/(2*a-b-c)**2+edjk**2 ) c c ea**2*DIF(SDK,a)**2+eb**2*DIF(SDK,b)**2+ec**2*DIF(SDK,c)**2+eddk** c $2*DIF(SDK,ddk)**2+edk**2*DIF(SDK,dk)**2 c esdk= dsqrt( $25*ddk**2*ea**2*(b-c)**2/(2*a-b-c)**4+25*ddk**2*eb**2*(a-c)**2/(2* $a-b-c)**4+25*ddk**2*ec**2*(a-b)**2/(2*a-b-c)**4+25*eddk**2*(b-c)** $2/(4*(2*a-b-c)**2)+edk**2 ) c ed1=eddj c c ea**2*DIF(D2,a)**2+eb**2*DIF(D2,b)**2+ec**2*DIF(D2,c)**2+eddk**2*D c $IF(D2,ddk)**2 c ed2=dsqrt( $ ddk**2*ea**2*(b-c)**2/(4.d0*(2.d0*a-b-c)**4) $ +ddk**2*eb**2*(a-c)**2/(4*(2*a-b-c)**4)+ddk**2*ec**2*(a-b)**2 $ /(4.d0*(2.d0*a-b-c)**4)+eddk**2*(b-c)**2 $ /(16.d0*(2.d0*a-b-c)**2) ) c c...errors in sextic constants c S=(2.D0*A-B-C)/(B-C) c c S:=(2*a-b-c)/(b-c) c c ea**2*DIF(S,a)**2+eb**2*DIF(S,b)**2+ec**2*DIF(S,c)**2 c es=dsqrt( $ 4.d0*(a**2*(eb**2+ec**2)-2.d0*a*(b*ec**2+c*eb**2)+b**2 $ *(ea**2+ec**2)-2.d0*b*c*ea**2+c**2*(ea**2+eb**2))/(b-c)**4 ) c c H3:=(hhk+16*d2*(sdk+2*(s**2-2)*d2)/(b-c))/(32*s**2/3+9) c c es**2*DIF(H3,s)**2+eb**2*DIF(H3,b)**2+ec**2*DIF(H3,c)**2+esdk**2*D c $IF(H3,sdk)**2+ehhk**2*DIF(H3,hhk)**2+ed2**2*DIF(H3,d2)**2 c eh3=dsqrt( $ 36864.d0*es**2*s**2*(b*hhk-c*hhk-d2*(91.d0*d2-16.d0*sdk))**2 $ /((b-c)**2*(32.d0*s**2+27.d0)**4)+2304.d0*d2**2*eb**2 $ *(2.d0*d2*(s**2-2)+sdk)**2/((b-c)**4*(32.d0*s**2+27.d0)**2) $ +2304.d0*d2**2*ec**2*(2.d0*d2*(s**2-2.d0)+sdk)**2/((b-c)**4 $ *(32.d0*s**2+27.d0)**2)+2304.d0*d2**2*esdk**2/((b-c)**2 $ *(32.d0*s**2+27.d0)**2)+9.d0*ehhk**2/(32.d0*s**2+27.d0)**2 $ +2304.d0*ed2**2*(4.d0*d2*(s**2-2.d0)+sdk)**2/((b-c)**2* $ (32.d0*s**2+27.d0)**2) ) c c H2:=(hhjk+10*h3+8*d2*(sdjk+2*s*d1+4*d2)/(b-c))/(4*s) c c ehhjk**2*DIF(H2,hhjk)**2+eh3**2*DIF(H2,h3)**2+ed2**2*DIF(H2,d2)**2 c $+esdjk**2*DIF(H2,sdjk)**2+ed1**2*DIF(H2,d1)**2+eb**2*DIF(H2,b)**2+ c $ec**2*DIF(H2,c)**2+es**2*DIF(H2,s)**2 c eh2= dsqrt( $ ehhjk**2/(16.d0*s**2)+25.d0*eh3**2/(4.d0*s**2) $ +4.d0*ed2**2*(2.d0*d1*s+8.d0*d2+sdjk)**2/(s**2*(b-c)**2) $ +4.d0*d2**2*esdjk**2/(s**2*(b-c)**2)+16.d0*d2**2*ed1**2 $ /(b-c)**2+4.d0*d2**2*eb**2*(2.d0*d1*s+4.d0*d2+sdjk)**2 $ /(s**2*(b-c)**4)+4.d0*d2**2*ec**2*(2.d0*d1*s+4.d0*d2+sdjk)**2 $ /(s**2*(b-c)**4)+es**2*(b*(10.d0*h3+hhjk)+c*(-10.d0*h3-hhjk) $ +8.d0*d2*(4.d0*d2+sdjk))**2/(16.d0*s**4*(b-c)**2) ) c c H1:=hhj-h3 c c ehhj**2*DIF(H1,hhj)**2+eh3**2*DIF(H1,h3)**2 c eh1= dsqrt( eh3**2+ehhj**2 ) c c SHJ:=hj-2*h2 c c ehj**2*DIF(SHJ,hj)**2+eh2**2*DIF(SHJ,h2)**2 c eshj= dsqrt( 4.d0*eh2**2+ehj**2 ) c c SHJK:=hjk+12*h2-16*s*h3 c c $ehjk**2*DIF(SHJK,hjk)**2+eh2**2*DIF(SHJK,h2)**2+es**2*DIF(SHJK,s)* c $*2+eh3**2*DIF(SHJK,h3)**2 ) c eshjk= dsqrt( 144.d0*eh2**2+256.d0*eh3**2*s**2+ehjk**2 $ +256.d0*es**2*h3**2 ) c c SHKJ:=hkj-10*h2+160*s/3*h3 c c ehkj**2*DIF(SHKJ,hkj)**2+eh2**2*DIF(SHKJ,h2)**2+es**2*DIF(SHKJ,s)* c $*2+eh3**2*DIF(SHKJ,h3)**2 c eshkj= dsqrt( $ 100.d0*eh2**2+25600.d0*eh3**2*s**2/9.d0+ehkj**2 $ +25600.d0*es**2*h3**2/9.d0) c c SHK:=hk-112*s/3*h3 c c ehk**2*DIF(SHK,hk)**2+es**2*DIF(SHK,s)**2+eh3**2*DIF(SHK,h3)**2 c eshk= dsqrt( $ 12544.d0*eh3**2*s**2/9.d0+ehk**2+12544.d0*es**2*h3**2/9.d0 ) C RETURN END C C_____________________________________________________________________________ C SUBROUTINE ERRSTA C C ERRors on S TO A conversion: C C EVALUATION OF ERRORS IN CALCULATED A-REDUCTION CONSTANTS FROM PROPAGATION C OF ERRORS ON THE INGOING S-REDUCTION CONSTANTS C C formulae from DERIVE and the generating relations for that program have C been left in the comment lines C IMPLICIT REAL*8 (A-H,O-Z) COMMON /DAT/AA(15),AS(15),AERROR(15),SERROR(15) C EQUIVALENCE (AA(1),A),(AA(2),B),(AA(3),C),(AA(4),DJ), * (AA(5),DJK),(AA(6),DK),(AA(7),DDJ),(AA(8),DDK), * (AA(9),HJ),(AA(10),HJK),(AA(11),HKJ),(AA(12),HK), * (AA(13),HHJ),(AA(14),HHJK),(AA(15),HHK) EQUIVALENCE (AS(1),SA),(AS(2),SB),(AS(3),SC),(AS(4),SDJ), * (AS(5),SDJK),(AS(6),SDK),(AS(7),D1),(AS(8),D2), * (AS(9),SHJ),(AS(10),SHJK),(AS(11),SHKJ),(AS(12),SHK), * (AS(13),H1),(AS(14),H2),(AS(15),H3) EQUIVALENCE (AERROR(1),EA),(AERROR(2),EB),(AERROR(3),EC), * (AERROR(4),EDJ),(AERROR(5),EDJK),(AERROR(6),EDK), * (AERROR(7),EDDJ),(AERROR(8),EDDK),(AERROR(9),EHJ), * (AERROR(10),EHJK),(AERROR(11),EHKJ),(AERROR(12),EHK), * (AERROR(13),EHHJ),(AERROR(14),EHHJK),(AERROR(15),EHHK) EQUIVALENCE (SERROR(1),ESA),(SERROR(2),ESB),(SERROR(3),ESC), * (SERROR(4),ESDJ),(SERROR(5),ESDJK),(SERROR(6),ESDK), * (SERROR(7),ED1),(SERROR(8),ED2),(SERROR(9),ESHJ), * (SERROR(10),ESHJK),(SERROR(11),ESHKJ),(SERROR(12),ESHK), * (SERROR(13),EH1),(SERROR(14),EH2),(SERROR(15),EH3) C c c...errors in quadratics c c c SIGMA:=(2*sa-sb-sc)/(sb-sc) c c A:=sa+10*d2 c B:=sb-4*(2*SIGMA+1)*d2 c C:=sc+4*(2*SIGMA-1)*d2 c c esa**2*DIF(A,sa)**2+esb**2*DIF(A,sb)**2+esc**2*DIF(A,sc)**2+ed2**2 c $*DIF(A,d2)**2 c c esa**2*DIF(B,sa)**2+esb**2*DIF(B,sb)**2+esc**2*DIF(B,sc)**2+ed2**2 c $*DIF(B,d2)**2 c c esa**2*DIF(C,sa)**2+esb**2*DIF(C,sb)**2+esc**2*DIF(C,sc)**2+ed2**2 c $*DIF(C,d2)**2 c ea=dsqrt( 100.d0*ed2**2+esa**2 ) c eb=dsqrt( $ (256.d0*d2**2*(esa**2*(sb-sc)**2+esb**2*(sa-sc)**2+ $ esc**2*(sa-sb)**2)+32.d0*d2*esb**2*(sb-sc)**2*(sa-sc)+ $ (sb-sc)**2*(16.d0*ed2**2*(4.d0*sa-sb-3.d0*sc)**2+esb**2 $ *(sb-sc)**2))/(sb-sc)**4 ) c ec=dsqrt( $ (256.d0*d2**2*(esa**2*(sb-sc)**2+esb**2*(sa-sc)**2+ $ esc**2*(sa-sb)**2)+32.d0*d2*esc**2*(sb-sc)**2*(sa-sb)+ $ (sb-sc)**2*(16.d0*ed2**2*(4.d0*sa-3.d0*sb-sc)**2+ $ esc**2*(sb-sc)**2))/(sb-sc)**4 ) c c...Quartics c c DJ:=sdj-2*d2 c DJK:=sdjk+12*d2 c DK:=sdk-10*d2 c DDJ:=(-d1) c DDK:=(-4*d2*SIGMA) c edj=dsqrt( 4.d0*ed2**2+esdj**2 ) edjk=dsqrt( 144.d0*ed2**2+esdjk**2 ) edk=dsqrt( 100.d0*ed2**2+esdk**2 ) eddj=ed1 eddk=dsqrt( $ d2**2*(64.d0*esa**2/(sb-sc)**2+64.d0*esb**2*(sa-sc)**2 $ /(sb-sc)**4+64.d0*esc**2*(sa-sb)**2/(sb-sc)**4) $ +16.d0*ed2**2*(2.d0*sa-sb-sc)**2/(sb-sc)**2 ) c c...Diagonal sextics c c HJ:=shj+2*h2 C C eshj**2*DIF(HJ,shj)**2+eh2**2*DIF(HJ,h2)**2 C ehj=dsqrt( 4.d0*eh2**2+eshj**2 ) C C HJK:=shjk-12*h2+16*SIGMA*h3 C c eshjk**2*DIF(HJK,shjk)**2+eh2**2*DIF(HJK,h2)**2+eh3**2*DIF(HJK,h3) C $**2+esa**2*DIF(HJK,sa)**2+esb**2*DIF(HJK,sb)**2+esc**2*DIF(HJK,sc) C $**2 C ehjk=dsqrt ( $ 144.d0*eh2**2+256.d0*eh3**2*(2.d0*sa-sb-sc)**2/(sb-sc)**2 $ +1024.d0*esa**2*h3**2/(sb-sc)**2+1024.d0*esb**2*h3**2*(sa-sc)**2 $ /(sb-sc)**4+1024.d0*esc**2*h3**2*(sa-sb)**2/(sb-sc)**4+eshjk**2 ) C c HKJ:=shkj+10*h2-160*SIGMA/3*h3 c c eshkj**2*DIF(HKJ,shkj)**2+eh2**2*DIF(HKJ,h2)**2+eh3**2*DIF(HKJ,h3) c $**2+esa**2*DIF(HKJ,sa)**2+esb**2*DIF(HKJ,sb)**2+esc**2*DIF(HKJ,sc) c $**2 c ehkj=dsqrt( $ 100.d0*eh2**2+25600.d0*eh3**2*(2.d0*sa-sb-sc)**2 $ /(9.d0*(sb-sc)**2) $ +102400.d0*esa**2*h3**2/(9.d0*(sb-sc)**2)+102400.d0*esb**2*h3**2 $ *(sa-sc)**2/(9.d0*(sb-sc)**4)+102400.d0*esc**2*h3**2*(sa-sb)**2 $ /(9.d0*(sb-sc)**4)+eshkj**2 ) c c HK:=shk+112*SIGMA/3*h3 c c eshk**2*DIF(HK,shk)**2+eh3**2*DIF(HK,h3)**2+esa**2*DIF(HK,sa)**2+e c $sb**2*DIF(HK,sb)**2+esc**2*DIF(HK,sc)**2 c ehk=dsqrt( $ 12544.d0*eh3**2*(2.d0*sa-sb-sc)**2/(9.d0*(sb-sc)**2) $ +50176.d0*esa**2*h3**2/(9.d0*(sb-sc)**2)+50176.d0*esb**2*h3**2 $ *(sa-sc)**2/(9.d0*(sb-sc)**4)+50176.d0*esc**2*h3**2*(sa-sb)**2/ $ (9.d0*(sb-sc)**4)+eshk**2 ) c c...off-diagonal sextics c c HHJ:=h1+h3 c c eh1**2*DIF(HHJ,h1)**2+eh3**2*DIF(HHJ,h3)**2 c ehhj=dsqrt( eh1**2+eh3**2 ) c c HHJK:=4*SIGMA*h2-10*h3-8*d2*(sdjk+2*SIGMA*d1+4*d2)/(sb-sc) c c eh2**2*DIF(HHJK,h2)**2+ c eh3**2*DIF(HHJK,h3)**2+ C ed2**2*DIF(HHJK,d2)**2+ c ed1**2*DIF(HHJK,d1)**2+ c esdjk**2*DIF(HHJK,sdjk)**2+ c esa**2*DIF(HHJK,sa)**2+ c esb**2*DIF(HHJK,sb)**2+ c esc**2*DIF(HHJK,sc)**2 c ehhjk= dsqrt( $ 16.d0*eh2**2*(2.d0*sa-sb-sc)**2/(sb-sc)**2 + $ 100.d0*eh3**2 + 64.d0*ed2**2*(2.d0*d1*(2.d0*sa-sb-sc) $ +(sb-sc)*(8.d0*d2+sdjk))**2/(sb-sc)**4 + $ 256.d0*d2**2*ed1**2*(2.d0*sa-sb-sc)**2/(sb-sc)**4 + $ 64.d0*d2**2*esdjk**2/(sb-sc)**2 + $ 64.d0*esa**2*(4.d0*d1*d2+h2*(sc-sb))**2/(sb-sc)**4 + $ 64.d0*esb**2*(2.d0*d1*d2*(4.d0*sa-sb-3.d0*sc)+(sb-sc) $ *(4.d0*d2**2+d2*sdjk+h2*(sc-sa)))**2/(sb-sc)**6 + $ 64.d0*esc**2*(2.d0*d1*d2*(4.d0*sa-3.d0*sb-sc)+(sb-sc) $ *(4.d0*d2**2+d2*sdjk+h2*(sb-sa)))**2/(sb-sc)**6 ) c c HHK:=(32*SIGMA**2/3+9)*h3-16*d2*(sdjk+2*(SIGMA**2-2)*d2)/(sb-sc) c c eh3**2*DIF(HHK,h3)**2+ed2**2*DIF(HHK,d2)**2+esa**2*DIF(HHK,sa)**2+ c $esb**2*DIF(HHK,sb)**2+esc**2*DIF(HHK,sc)**2+esdjk**2*DIF(HHK,sdjk) c $**2 c ehhk=dsqrt( $ (9216.d0*d2**4*(16.d0*esa**2*(sb-sc)**2*(2.d0*sa-sb-sc)**2+esb**2 $*(12.d0*sa**2-8.d0*sa*(sb+2.d0*sc)-sb**2+sc*(10.d0*sb+3.d0*sc))**2 $ +esc**2*(12.d0*sa**2-8.d0*sa*(2.d0*sb+sc)+3.d0*sb**2+10.d0*sb*sc $ -sc**2)**2)+9216.d0*d2**3*sdjk*(sb-sc)**2*(esb**2*(12.d0*sa**2 $ -8.d0*sa*(sb+2.d0*sc)-sb**2+sc*(10.d0*sb+3.d0*sc))+esc**2*(12.d0 $ *sa**2-8.d0*sa*(2.d0*sb+sc)+3.d0*sb**2+10.d0*sb*sc-sc**2))+768.d0 $ *d2**2*(sb-sc)*(48.d0*ed2**2*(sb-sc)*(4.d0*sa**2-4.d0*sa*(sb+sc) $ -sb**2+sc*(6.d0*sb-sc))**2-128.d0*esa**2*h3*(sb-sc)**2*(2.d0*sa $ -sb-sc)**2+esb**2*(3.d0*sdjk**2*(sb-sc)**3-32.d0*h3*(2.d0*sa**2 $ +sa*(-sb-3.d0*sc)+sc*(sb+sc))*(12.d0*sa**2-8.d0*sa*(sb+2.d0*sc) $ -sb**2+sc*(10.d0*sb+3.d0*sc)))+esc**2*(3.d0*sdjk**2*(sb-sc)**3 $ -32.d0*h3*(2.d0*sa**2+sa*(-3.d0*sb-sc)+sb*(sb+sc))*(12.d0*sa**2 $ -8.d0*sa*(2.d0*sb+sc)+3.d0*sb**2+10.d0*sb*sc-sc**2))+3.d0 $ *esdjk**2*(sb-sc)**5)+6144.d0*d2*sdjk*(sb-sc)**3*(3.d0*ed2**2 $ *(sb-sc)*(4.d0*sa**2-4.d0*sa*(sb+sc)-sb**2+sc*(6.d0*sb-sc))-2.d0 $ *h3*(2.d0*sa-sb-sc)*(esb**2*(sa-sc)+esc**2*(sa-sb)))+(sb-sc)**2 $ *(2304.d0*ed2**2*sdjk**2*(sb-sc)**4+eh3**2*(sb-sc)**2*(128.d0 $ *sa**2-128.d0*sa*(sb+sc)+59.d0*sb**2+10.d0*sb*sc+59.d0*sc**2)**2 $ +16384.d0*h3**2*(2.d0*sa-sb-sc)**2*(esa**2*(sb-sc)**2+esb**2 $ *(sa-sc)**2+esc**2*(sa-sb)**2)))/(9.d0*(sb-sc)**8) ) C RETURN END C C_____________________________________________________________________________ c SUBROUTINE CONFOR(CONVAL,ERVAL,NDCON,NDEROR,NERD) C C Parameter and error formatting for output c c CONVAL - String containing the parameter value. This is to be c input in F format and will be replaced on output by the result c string of length extending to the last digit of the error c ERVAL - String containing the error value. This is to be c input in F format and will be replaced on output by the result c string, which does not contain the decimal point and is meant c to be included in brackets c NDCON - The number of digits in the CONVAL string (inclusive of any c leading zeros) c NDEROR - The number of digits in the ERVAL string (just the significant c digits) and is either equal to NERD, or is larger if there c are more significant digits than NERD before the decimal point. c NERD - the number of desired error digits, which is to be set on input. c c Both NDCON and NDEROR are generated on output C CHARACTER(27) CONVAL,ERVAL,CONOUT,EROUT C NDEROR=0 NDNOTZ=0 ICZERO=0 erout='???????????' C C...Go through digits of parameter value adding those to output buffer while C checking at the same time digits of the error value, reacting as C necessary. C Terminate when either NERD digits in the error value are C transferred or, if error has more digits before the decimal C point, the decimal point is reached. C conout(1:2)=conval(1:2) ndcon=2 DO 1 N=3,27 NDCON=NDCON+1 CONOUT(NDCON:NDCON)=CONVAL(N:N) C C...ensure that zero precedes the decimal point and use ICZERO to monitor C whether decimal point has been reached IF(CONVAL(N:N).EQ.'.')ICZERO=1 IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.' ') * CONOUT(NDCON-1:NDCON-1)='0' IF(CONVAL(N:N).EQ.'.'.AND.CONVAL(N-1:N-1).EQ.'-')THEN CONOUT(NDCON-2:NDCON-2)='-' CONOUT(NDCON-1:NDCON-1)='0' ENDIF C C...use NDNOTZ (number of digits not zero) to monitor whether significant C digits in parameter value have been reached IF(CONVAL(N:N).GE.'1'.AND.CONVAL(N:N).LE.'9') * NDNOTZ=NDNOTZ+1 C C...do not transfer error digit if it is a leading zero, dot or space IF(NDEROR.EQ.0 .AND. (ERVAL(N:N).EQ.' '.OR. * ERVAL(N:N).EQ.'0'.OR.ERVAL(N:N).EQ.'.') )GOTO 1 C C...exit if error larger than value and decimal point reached IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 C C...do not transfer the dot in error value IF(ERVAL(N:N).EQ.'.')GOTO 1 C C...transfer error digits until NERD or, if more, the first significant C digit in value is reached NDEROR=NDEROR+1 EROUT(NDEROR:NDEROR)=ERVAL(N:N) IF(NDEROR.GE.NERD .AND. NDNOTZ.GT.0 .AND. ICZERO.EQ.1)GOTO 2 1 CONTINUE c c...output string with fitted value (rounding is to be carried out externally) c 2 CONVAL(1:NDCON)=CONOUT(1:NDCON) c c...output string with error (rounded if necessary) c if(nderor.gt.0)then c write(*,'(1x,2a,5x,a)') DEBUG c * 'ERROR: ',erval,erout(1:nderor) DEBUG c if(erval(ndcon+1:ndcon+1).eq.'5'.or. * erval(ndcon+1:ndcon+1).eq.'6'.or. * erval(ndcon+1:ndcon+1).eq.'7'.or. * erval(ndcon+1:ndcon+1).eq.'8'.or. * erval(ndcon+1:ndcon+1).eq.'9')then read(erout(1:nderor),*)ieror net=int(dlog10(dble(ieror))) ! Fortran2018 ieror=ieror+1 net1=int(dlog10(dble(ieror))) ! Fortran2018 if(net1-net.gt.0)nderor=nderor+1 write(erout,'(i12)')ieror erval(1:nderor)=erout(12-nderor+1:12) else ERVAL(1:NDEROR)=EROUT(1:NDEROR) endif endif C RETURN END C C------------------------------------------------------------------------- C-------------------------------------------------------------------------