PROGRAM NEWMEOH ************************************************************************ * * * Deconvolution of the measured absorbances with a kinetic mechanism * * specified in subroutine KINETICS for the solvation of electrons. * * (The actual model is that from Rad. Phys. Chem. ............... * * The pulse shape used is taken as numerical data from "pulse.dat". * * (The norm of the pulse (overall sum) is 100.00000) * * * * Evaluation of MONOCANAL METHANOL data at 13 different WAVELENGTHs * * * * The spectral shape (absorbances of the 4 species at each wave- * * length) is contained in the file "spect.txt". Relative amplitude * * of the absorbance spectra is specified by the parameters K1...K4 * * in the file "mall.new", along with the characteristic times. * * (In this model, the 5th species' (solvated electron) spectrum is * * taken from independent experiment and fixed in SUBROUTINE MODEL.) * * * * To be linked with the MMRQT routine. * * * ************************************************************************ CCCCC Modified to run on RISC/6000 workstations CCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IMPLICIT DOUBLE PRECISION (D) REAL*8 P0,PP,W COMMON G(1002),DABS(13,4) DIMENSION DT(13,300,1),DV(13,300,1),DN(13,300,1),P0(15),PP(15), 1 IP(15),NM(13),NX(13),NY(13),W(9,9),A(10) CHARACTER NAME*40,FNAME*20,FFIT*20,FPRT*20,SHIT*20 * FPRT='prt.dat' FFIT='fit.dat' * ******************************************************* OPEN(7,FILE='mall.new',STATUS='OLD') ******************************************************* READ(7,100)NAME WRITE(*,100)NAME READ(7,100)SHT WRITE(*,100)SHT READ(7,100)SHT WRITE(*,100)SHT 100 FORMAT(A40) 150 FORMAT(1X,/,/,1X,A40,/) * READ(7,*)NN,NMM,NXM,NYM,NS,MW,MD,NSTEP,DSTEP WRITE(*,*)NN,NMM,NXM,NYM,NS,MW,MD,NSTEP,DSTEP * READ(7,100)SHT READ(7,*)(PP(I),I=1,NS) READ(7,*)(IP(I),I=1,NS) WRITE(*,100)SHT WRITE(*,*)(PP(I),I=1,NS) WRITE(*,*)(IP(I),I=1,NS) NP=0 DO 990 I=1,NS IF (IP(I).EQ.1) THEN NP=NP+1 P0(NP)=PP(I) IP(NP)=I ENDIF 990 CONTINUE * 230 FORMAT(1X,A40) * DO 500 K=1,NN READ(7,100)SHT READ(7,*)NM(K),NX(K),NY(K) WRITE(*,100)SHT WRITE(*,201)NM(K),NX(K),NY(K) 201 FORMAT(3I5) DO 500 I=1,NM(K) IF (MW.NE.5) GOTO 399 READ(7,*)(DT(K,I,J),J=1,NX(K)),(DV(K,I,J),J=1,NY(K)), 1 (DN(K,I,J),J=1,NY(K)) GOTO 500 399 READ(7,*)(DT(K,I,J),J=1,NX(K)),(DV(K,I,J),J=1,NY(K)) * WRITE(*,*)(DT(K,I,J),J=1,NX(K)),(DV(K,I,J),J=1,NY(K)) * 500 CONTINUE * DO 3988 J=1,NY(1) DO 3988 I=1,NM(1) 3988 DV(1,I,J)=7.766/7.528*DV(1,I,J) DO 3989 J=1,NY(2) DO 3989 I=1,NM(2) 3989 DV(2,I,J)=9.084/8.747*DV(2,I,J) DO 3990 J=1,NY(3) DO 3990 I=1,NM(3) 3990 DV(3,I,J)=10.369/10.250*DV(3,I,J) DO 3991 J=1,NY(11) DO 3991 I=1,NM(11) 3991 DV(11,I,J)=1.2*DV(11,I,J) DO 3992 J=1,NY(12) DO 3992 I=1,NM(12) 3992 DV(12,I,J)=5.0*DV(12,I,J) DO 3993 J=1,NY(13) DO 3993 I=1,NM(13) 3993 DV(13,I,J)=7.0*DV(13,I,J) * OPEN(3,FILE='pulse.dat',STATUS='OLD') READ(3,100)SHIT READ(3,100)SHIT READ(3,*)NSL READ(3,100)SHIT READ(3,610)(G(I),I=1,NSL) OPEN(4,FILE='spect.txt',STATUS='OLD') DO 1313 I=1,13 1313 READ(4,*)(DABS(I,J),J=1,4) 605 FORMAT(1H1) 610 FORMAT(8(F9.6,1X)) * * CALL MMRQT(NAME,NN,NMM,NXM,NYM,NM,NX,NY,NP,NS,MW,MD,P0,PP,IP, 1 DT,DV,DN,NSTEP,DSTEP,W,A,NSL,FFIT,FPRT) * WRITE(*,655) * *********************** model explanation ************************** * 655 FORMAT(//,5X, /, 1 5X,' Parameters: k0..k5, (a1.. a5) ') * > < *********************************************************************** * 300 CONTINUE WRITE(*,605) END * SUBROUTINE USERW(NN,MN,NMM,NXM,NYM,M,NX,NY,NP,NS,DT,DV,DW,DWS, 1 DP,AC,I) ************************************************************************ * Calculates USER DEFINED weights for the dependent variables. * * Input variables: NX= number of independent variables * * NY= number of dependent variables * * M= index of the actual measured point * * DV(M,I)= value of the ith dependent variable * * DW(I,I)= actual weight of the ith dep.var. * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION P,W DIMENSION DT(13,300,1),DV(13,300,1),DW(9,9),DP(15),AC(10) C *** Here you can write your own weight-function.************** * IF(M.LE.14) DW(1,1)=DWS/.3 IF(M.GT.14) DW(1,1)=DWS/1.5 * RETURN END SUBROUTINE MODEL(MN,NX,NY,NP,NS,X,Y,PT,PP,IP,B,NSL) IMPLICIT DOUBLE PRECISION (A,C,D,K-L,T,Z) DOUBLE PRECISION X,Y,P,PT,PP,PKIN,PKK14,NRCB,NRCB1 COMMON G(1002),DABS(13,4) DIMENSION X(15),Y(15),P(15),PP(NS),IP(NS),PT(15),B(10),C(10), 1 PKIN(15) * SL = REAL(NSL) DO 980 I=1,NS 980 P(I)=PP(I) DO 970 I=1,NP 970 P(IP(I))=PT(I) IF (MN.EQ.1) THEN A5 = 7.766 THALF = 0.410 T0 = 0.131 ENDIF IF (MN.EQ.2) THEN A5 = 9.084 THALF = 0.410 T0 =-0.193 ENDIF IF (MN.EQ.3) THEN A5 = 10.369 THALF = 0.675 T0 = 0.000005 ENDIF IF (MN.EQ.4) THEN A5 = 10.322 THALF = 0.675 T0 = 0.278 ENDIF IF (MN.EQ.5) THEN A5 = 9.623 THALF = 0.480 T0 =-0.186 ENDIF IF (MN.EQ.6) THEN A5 = 8.294 THALF = 0.480 T0 = 1.387 ENDIF IF (MN.EQ.7) THEN A5 = 6.560 THALF = 0.480 T0 = 0.026 ENDIF IF (MN.EQ.8) THEN A5 = 4.866 THALF = 0.480 T0 = 0.658 ENDIF IF (MN.EQ.9) THEN A5 = 3.465 THALF = 0.480 T0 = 0.170 ENDIF IF (MN.EQ.10) THEN A5 = 2.409 THALF = 0.480 T0 = 0.536 ENDIF IF (MN.EQ.11) THEN A5 = 0.4760 THALF = 0.480 T0 = 0.256 ENDIF IF (MN.EQ.12) THEN A5 = 0.000 THALF = 0.480 T0 = 0.098 ENDIF IF (MN.EQ.13) THEN A5 = 0.00 THALF = 0.480 T0 = 0.021 ENDIF DIV = 2. * THALF / SL TK = X(1) - T0 - THALF + DIV/2. Y(1) = 0.0 **** specific absorbances ********** * A1 = P(6) * DABS(MN,1) A2 = P(7) * DABS(MN,2) A3 = P(8) * DABS(MN,3) A4 = P(9) * DABS(MN,4) * * * k0: PKIN(1) = 1./ P(1) * k1: PKIN(2) = 1./ P(2) * k2: PKIN(3) = 1./ P(3) * k3: PKIN(4) = 1./ P(4) * k: PKIN(5) = 1./ P(5) * number of states: PKIN(6) = 19. * PKK14=PKIN(1)+PKIN(4) * ************************************ * **** Computation of the absorbances ********************************* * * DO 2222 NK=1,NSL T = TK + NK * DIV * * * IF(T.LT.0.)GOTO 2222 CALL KINETICS(PKIN,T,C) IF ((MMMN.EQ.19).OR.(MMMN.EQ.29)) THEN Y(1)=Y(1)+G(NK)*(A1*C(1)+A2*C(2)+A3*C(3) X +(A4*C(4)+A5*C(5)) X *NRCB1(T,1./PKIN(1))*PKIN(1)/PKK14 X +A5*C(5) X *NRCB1(T,1./PKIN(4))*PKIN(4)/PKK14)/100. ELSE Y(1)=Y(1)+G(NK)*(A1*C(1)+A2*C(2)+A3*C(3) X +(A4*C(4)+A5*C(5)) X *NRCB(T,1./PKIN(1))*PKIN(1)/PKK14 X +A5*C(5) X *NRCB(T,1./PKIN(4))*PKIN(4)/PKK14)/100. ENDIF 2222 CONTINUE * RETURN END * * SUBROUTINE KINETICS (PKIN,TIME,C) ************************************************************************ * * * Calculates the CONCENTRATIONS needed for the absorbance evaluation * * in SUBROUTINE MODEL. In that latter program, there is only calcul- * * ation of the absorbances at different steps of convolution for * * different wavelengths and associated pulse characteristics. * * TIME = the time parameter at the actual pulse * * C(I) = the ith concentration function at time TIME * * PKIN(15) contains the necessary kinetic parameters * * * ************************************************************************ IMPLICIT DOUBLE PRECISION (A,C,D,K-L,S,T,Z) DOUBLE PRECISION PKIN DIMENSION C(10),PKIN(15),KF(20),KSF(20),K1F(20),K2F(20), 1 KH1(20),KH2(20) * C(1) = 0.0D0 C(2) = 0.0D0 C(3) = 0.0D0 C(4) = 0.0D0 C(5) = 0.0D0 IF(TIME.LT.0) GOTO 666 ******* PART I. Assignement of the kinetic parameters *********** * K0 = PKIN(1) K1 = PKIN(2) K2 = PKIN(3) K3 = PKIN(4) K = PKIN(5) NST = PKIN(6) * KS = K0 + K3 KKS = K - KS CCC IF(ABS(KKS).LE.1.0D-30) KKS = SIGN(1.0D-30,KKS) K1S = K1 - KS CCC IF(ABS(K1S).LE.1.0D-30) K1S = SIGN(1.0D-30,K1S) K2S = K2 - KS CCC IF(ABS(K2S).LE.1.0D-30) K2S = SIGN(1.0D-30,K2S) KK1 = K - K1 CCC IF(ABS(KK1).LE.1.0D-30) KK1 = SIGN(1.0D-30,KK1) KK2 = K - K2 CCC IF(ABS(KK2).LE.1.0D-30) KK1 = SIGN(1.0D-30,KK2) K21 = K2 - K1 CCC IF(ABS(K21).LE.1.0D-30) K21 = SIGN(1.0D-30,K21) * KN = (K/KKS)**DBLE(NST) T = TIME KT = K*T KST = KS*T K1T = K1*T K2T = K2*T KKST = KKS*T KK1T = KK1*T KK2T = KK2*T * CCC IF(KT.GE.180.D0) THEN CCC KE = 0.D0 CCC ELSE KE = DEXP(-KT) CCC ENDIF * CCC IF(KST.GE.180.D0) THEN CCC KSE = 0.D0 CCC ELSE KSE = DEXP(-KST) CCC ENDIF * CCC IF(K1T.GE.180.D0) THEN CCC K1E = 0.D0 CCC ELSE K1E = DEXP(-K1T) CCC ENDIF * CCC IF(K2T.GE.180.D0) THEN CCC K2E = 0.D0 CCC ELSE K2E = DEXP(-K2T) CCC ENDIF * KDE1 = (K1E - K2E) / K21 KDE2 = (KSE - K2E) / K2S KF(1) = 1 KSF(1) = 1 K1F(1) = 1 K2F(1) = 1 KH1(1) = 1 / KK1 KH2(1) = 1 / KK2 DO 131 I = 2,NST KF(I) = KF(I-1) * KT / (I-1) KSF(I) = KSF(I-1) * KKST / (I-1) K1F(I) = K1F(I-1) * KK1T / (I-1) K2F(I) = K2F(I-1) * KK2T / (I-1) KH1(I) = KH1(I-1) * KKS / KK1 KH2(I) = KH2(I-1) * KK1 / KK2 131 CONTINUE C******** PART II. Computation of the concentration functions ********* * SUM1 = 0.D0 SUM2 = 0.D0 SUM3 = 0.D0 SUM4 = 0.D0 DO 1 I = 0,NST-1 SUMB3 = 0.D0 SUMB4 = 0.D0 DO 2 J = 0,I SUMBB4 = 0.D0 DO 3 M = 0,J SUMBB4 = SUMBB4 + K2F(M+1) 3 CONTINUE SUMB3 = SUMB3 + K1F(J+1) SUMB4 = SUMB4 + KH2(J+1) * (K2E - KE * SUMBB4) 2 CONTINUE SUM1 = SUM1 + KF(I+1) SUM2 = SUM2 + KSF(I+1) SUM3 = SUM3 + KH1(I+1) * (K1E - KE * SUMB3) SUM4 = SUM4 + KH1(I+1) * (KDE1 - SUMB4) 1 CONTINUE C(1) = KE * SUM1 C(2) = KN * (KSE - KE * SUM2) C(3) = KN * K0 * ( 1 / K1S * (KSE - K1E) - SUM3) C(4) = KN * K0 * K1 * ( 1 / K1S * (KDE2 - KDE1) - SUM4) C(5) = 1.D0 - C(1) - C(2) - C(3) - C(4) * ******** End concentration function computation *************** * 666 RETURN END * * * * FUNCTION NRCB(T,TLOC) ************************************************************************ * Calculates the NONRECOMBINED RATIO for excess electrons in METHANOL * * as a function of time elapsed from the localization of the electron * * after ionization. The function is a sum of a constant A0 and three * * exponential terms with preexponentials Ai and exponents Si so that * * there is a ratio RCMB(100 ps)/RCMB(40 ps) = 0.93, which is the ex- * * perimentally observed one. The function is an approximation of what * * has been found in simulations by Thomas Goulet. * * ( From above 540 nm ) * * The expected value of this exponential sum is then calculated over * * the exponential distribution of the localized electrons' concent- * * ration function. * ************************************************************************ DOUBLE PRECISION Q,R,NRCB,T,TLOC DIMENSION A(3),S(3),Q(3),R(3) * ** Constants of the exponential-sum fit to simulated recombination ** * A0 = .609 A(1) = .15171 A(2) = .23179 A(3) = .007 S(1) = 7.1704 S(2) = 169.47 S(3) = 31.055 * ** Computation of the localization-time-corrected function constants ** * DO 100 I = 1,3 Q(I) = A(I)*S(I)/(S(I)-TLOC) R(I) = (S(I)-TLOC)/(S(I)*TLOC) 100 CONTINUE * ***** Computation of the NON-RECOMBINED ratio at time T ****** * IF(T.LE.0.D0) THEN NRCB = 0. ELSE NRCB = A0 DO 200 J = 1,3 200 NRCB = NRCB + 1 Q(J)*(1.-DEXP(-R(J)*T))*DEXP(-T/S(J))/(1.-DEXP(-T/TLOC)) ENDIF * RETURN END * * FUNCTION NRCB1(T,TLOC) ************************************************************************ * Calculates the NONRECOMBINED RATIO for excess electrons in METHANOL * * as a function of time elapsed from the localization of the electron * * after ionization. The function is a sum of a constant A0 and three * * exponential terms with preexponentials Ai and exponents Si so that * * there is a ratio RCMB(100 ps)/RCMB(40 ps) = 0.92, which is the ex- * * perimentally observed one. The function is an approximation of what * * has been found in simulations by Thomas Goulet. * * ( From 500 to 540 nm ) * * The expected value of this exponential sum is then calculated over * * the exponential distribution of the localized electrons' concent- * * ration function. * ************************************************************************ DOUBLE PRECISION Q,R,NRCB1,T,TLOC DIMENSION A(3),S(3),Q(3),R(3) * ** Constants of the exponential-sum fit to simulated recombination ** * A0 = .56856 A(1) = .18257 A(2) = .24214 A(3) = .00673 S(1) = 6.971 S(2) = 153.474 S(3) = 17.25939 * ** Computation of the localization-time-corrected function constants ** * DO 100 I = 1,3 Q(I) = A(I)*S(I)/(S(I)-TLOC) R(I) = (S(I)-TLOC)/(S(I)*TLOC) 100 CONTINUE * ***** Computation of the NON-RECOMBINED ratio at time T ****** * IF(T.LE.0.D0) THEN NRCB1 = 0. ELSE NRCB1 = A0 DO 200 J = 1,3 200 NRCB1 = NRCB1 + 1 Q(J)*(1.-DEXP(-R(J)*T))*DEXP(-T/S(J))/(1.-DEXP(-T/TLOC)) ENDIF * RETURN END * ************************************************************************ * * SUBROUTINE DERJACOBI(NM,NX,NY,NP,NS,X,Y,DG,P,PP,IP,AC,IC) ************************************************************************ * Calculates the Jacobian matrix of the derivatives dYi/dPj. * * If you wish to use analytical derivatives, put flag DM = 1 in * * the driver program of MMRQT and write this subprogram to evaluate * * the Jacobian matrix DG(i,j) = dY(i)/dP(j). * ************************************************************************ IMPLICIT DOUBLE PRECISION (D) DOUBLE PRECISION X,Y,P DIMENSION X(NX),DG(15,15),P(NS),PP(NS),IP(NS),Y(NY),AC(10) C *** Write here the evaluation of DG ********************************* * * * * ************************************************************************ CCCCC write(*,*)(dg(1,j),j=1,4) RETURN END