$debug C============================================================================= C C Program for fitting rotational spectra in excited states of a doubly C degenerate vibrational mode in C3v symmetric top C C written by J.Demaison, original name diagc3v c c cosmetics/modifications by Z.Kisiel ver 6.III.1997 C------------------------------------------------------------------------------ c c files: ve.inp - default input c ve.out - default ouput c ve.hdr - new header for input file c c Note that an additional line is now at the top of the input c file which contains comment about the data set c C------------------------------------------------------------------------------ c constants: q = -q (LFITD) c q.j = -q.J (LFITD) c c For v=1, kl-1=0 the K,l assignments of low,high frequency components c alternates with J ie. assignment follows A1, A2 transition symmetry c c For explanation of the statistical tests see J.Demaison, J.Cosleou, c R.Bocquet, and A.G.Lesarri, J.Mol.Spectrosc. 167,400-418(1994) c c t(i) = "jacknifed residuals" c hii = diagonal elements of the H matrix, highest values identify the most c 'influential' points, ie those that predominantly determine some c constant c DFBETASj(i) = highvalues identify sensitivity of evaluation of parameter c j to dropping data point i c c C------------------------------------------------------------------------------ c c Switches and key variables: c c top line in data: v, itmax, iw where c v = integer being the number of quanta of excitation c itmax = number of iterations c iw = 0 apply equal weights to data c = 1 allow weighting according to errors in input c c icst = number of fixed parameters c numl = number of lines in the fit c c Use of arrays: c c c(i) - constants c ilib(i) - integer fitting parameter (1/0=floating/fixed) of ith constant c lib(i) - character fitting parameters ( or *) c icc(i) - pointers to ICST fixed parameters c jn(i) - J quantum number of line i c kn(i) - K quantum number of line i c ln(i) - l quantum number of line i c fr(i) - measured frequency of line i c w(i) - weight of line i, being an inverse of 0.001*(input error) c da(i) - at the end of the main iteration loop contains c corrections to constants resulting from the iteration such c that c.i+1=c.i-da c y(j) - standard deviation of j'th fitted constant c a(i,j), a1(i,j), a2(i,j), e1(i,j), t1(i,j) c - energy matrices of the Hamiltonian of size ca. c (v+1)(2J+1)/3 c amat(i,j) - derivative of line i with respect to constant j, some are c evaluated numerically, others analytically c C============================================================================= c c implicit none c parameter (maxlin=330,nconst=27,MAXDIM=100) c logical lg4,lg(0:75) ! jmax=75 character coment*79,testi*3 character*1 lib(nconst) character*6 param(nconst) c integer v,j,na,n1,ne1,nad,na1,na2,l,k,ll,l2,kk,k2,kl,tau,it, + ia,ie,ja,je,jj,indice,jmax,numl,iw,kle, + ka(maxdim),la(maxdim),ke(maxdim),le(maxdim),jn(maxlin), + kn(maxlin),ln(maxlin),indla(maxdim),indle(maxdim), + itm,icst,itmax,icc(nconst),mm,ilib(nconst),nq,nx c integer i,m,n,ib,ic,in,ip,ir,jm,nrot,counta,countb real*4 hat,hat0,hat1 c real*8 a(maxdim,maxdim),e1(maxdim,maxdim),t1(maxdim,maxdim), + a1(maxdim,maxdim),a2(maxdim,maxdim),c(nconst), + da(maxdim),va(maxdim,maxdim),de(maxdim),ve(maxdim,maxdim), + s,fr(maxlin),w(maxlin), + pp,f1,f2,en(maxlin,2),fber(maxlin),amat(maxlin,nconst), + deriva(maxdim,6),derive(maxdim,6),ai(maxdim,maxdim), + verb(maxlin),sdev(maxlin),y(nconst),h1,h2,ec,varz,dacc c real*4 rint(maxlin) real*8 wk(maxlin) integer*2 ipt(maxlin) common /DERIV/amat COMMON /sort/WK,ipt common /work/a,t1 c data param/' B ',' DJ ',' DJK ',' HJ ',' HJK ',' HKJ ', + ' nJ ',' nJJ ',' nJK ',' nJJJ ',' nJJK ',' nJKK ',' B3 ', + ' q¡ ',' qj ',' qjj ',' qk ',' Az ',' X ',' A ', + ' DK ',' nK ',' HK ',' B5J ',' B5K ',' B5l ',' emp '/ spnucl=1.5 c open(1,file='v2e.out',status='unknown') open(2,file='v2e.inp',status='old') c read(2,'(a)')coment read(2,*)v,itmax,iw write(6,2000)coment write(1,2000)coment write(6,*)' v =',v,' iter =',itmax,' weight =',iw write(1,*)' v =',v,' iter =',itmax,' weight =',iw 2000 format(1x,78(1h-)/1x,a/1x,78(1h-)/) itmax=itmax+1 c c Constants and their numbers: c c 1 2 3 4 5 6 7 8 9 10 11 12 c B, DJ, DJK, HJ, HJK,HKJ, nJ, nJJ, nJK, nJJJ, nJJK, nJKK, c c 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 c B3, q, q1, q2, qK, Az, X, A, DK, nK, HK B5J B5K B5l emp c c The 27th parameter is an empirical one which can be used for tests c icst=0 ! read constants, fixed if ilib(i)=0 do i=1,nconst read(2,*)c(i),ilib(i) if(ilib(i).eq.0) then icst=icst+1 icc(icst)=i lib(i)=' ' else lib(i)='*' end if ! ilib(i)=0 end do ! i c nq=2 nx=6 ! limit on numerical evaluation of derivatives if(ilib(17).eq.0)nq=1 if((ilib(14)+ilib(15)+ilib(16)+ilib(17)).eq.0)nq=0 if((ilib(19)+ilib(20)+ilib(22)).eq.0)nx=3 if((ilib(18)+ilib(19)+ilib(20)+ilib(22)).eq.0)nx=0 f2=0.d0 ! read transition frequencies do i=1,maxlin f1=1.d0 2006 read(2,*,err=2006)jn(i),kn(i),ln(i),fr(i),pp if(jn(i).eq.0)goto 2006 if(jn(i).eq.1111)exit if(iw.eq.1)f1=f1/(pp*1.d-3) ! weighting sdev(i)=0.d0 f2=f2+f1*f1 w(i)=f1 end do numl=i-1 hat0=float(nconst-icst)/numl f2=dsqrt(numl/f2) write(6,1005)numl,f2,hat0 write(1,1005)numl,f2,hat0 do i=1,numl w(i)=w(i)*f2 end do varz=0.d0 c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c c...Beginning of iteration loop, executed ITMAX times c do itm=1,itmax c write(6,2001)itm-1 write(1,2001)itm-1 2001 format(1x/' Iteration no:',i2,' --> Constant ', * ' Deviation Change'/) ncount=0 do i=1,nconst if(c(i).ne.0.d0)then if(lib(i).eq.'*')then ncount=ncount+1 write(6,1001)param(i),lib(i),c(i),y(ncount),da(i) write(1,1001)param(i),lib(i),c(i),y(ncount),da(i) else write(6,1001)param(i),lib(i),c(i) write(1,1001)param(i),lib(i),c(i) endif ! lib(i) endif ! c(i) end do ! i c lg4=itm.eq.itmax dacc=0.d0 hat1=hat0*varz indice=1 do i=1,75 ! initialisation lg(i)=.false. end do do i=1,numl fber(i)=0.d0 do j=1,nconst amat(i,j)=0.d0 end do ! j end do ! i c c...go over all lines, first work on J'', then J' of the line (selected by in) c do i=1,numl j=jn(i) if(j.lt.0) cycle do in=1,2 j=j+in-1 if(lg(j).eqv..false.) then lg(j)=.true. c c...calculate energy levels, array sizes are: c c na=rang de A c ne1=rang de E1 c na1=rang de A+ c na2=rang de A- c jmax=j c do while(jmax.gt.kmax) ! limitation de la taille des matrices c jmax=jmax-2 c end do na=(v+1)*(2*jmax+1) n1=mod(na,3) c select case(n1) case(0) na=na/3 ne1=na case(1) ne1=na/3 na=ne1+1 case(2) na=na/3 ne1=na+1 case default end select c n1=mod(na,2) nad=na/2 na1=nad na2=nad c c...check of legal array size c if(na.gt.maxdim.or.ne1.gt.maxdim.or.na1.gt.maxdim.or. * na2.gt.maxdim)then write(*,2007)maxdim,j,na,ne1,na1,na2 2007 format(1x//' ERROR: maximum allowed array size of',i4, * ' exceeded for:'//' J =',i4,' na,ne1,na1,na2: ', * 4(i5)//' Increase MAXDIM in PARAMETER statements?'//) stop endif c if(n1.eq.1) then if(mod(jmax,2).eq.0) then na1=na1+1 else na2=na2+1 end if end if c c...initialisation des matrices c do ia=1,na do ie=1,na a(ia,ie)=0.d0 t1(ia,ie)=0.d0 end do end do do ia=1,ne1 do ie=1,ne1 e1(ia,ie)=0.d0 end do end do c c...set up energy matrix for a given J, for all values of v as counted by l c ja=0 ia=0 ie=0 jj=j*(j+1) do l=-v,v,2 ll=l+2 l2=l*l do k=-jmax,jmax tau=mod(k-l,3) if(tau.lt.0)tau=tau+3 if(tau.eq.2) cycle ! matrice Eb k2=k*k kl=k*l kk=k+2 c select case(tau) c case(0) ! k-l=3n -> A+ & A- ia=ia+1 t1(ia,ia)=1.d0 if(na1.ge.na2) then ! stockage de k et l dans ka et la if(ia.le.na1) then ka(ia)=k la(ia)=l else if(((kl.eq.1).and.(v.eq.1)).or.((kl.eq.4)) + .and.(v.eq.2)) then indice=ia-na1 ! splitting A1-A2 ka(na1+1)=k la(na1+1)=l end if ! kl end if ! ia else if(ia.le.na1) then if(kl.eq.4) then ! splitting A1-A2 indice=ia ka(na2+1)=k la(na2+1)=l end if ! kl else ka(ia-na1)=k la(ia-na1)=l end if ! ia end if ! na1 c deriva(ia,3)=-2.d0*kl ! deriv Az deriva(ia,4)=l2 ! deriv X deriva(ia,5)=k2 ! deriv A deriva(ia,6)=k2*kl ! deriv eta.K a(ia,ia)=c(19)*l2+(c(1)+c(13)*l2-c(3)*k2)*jj+ @ (c(20)-c(1)-c(13)*l2+(c(23)*k2-c(21)+c(6)*jj)*k2)*k2+ @ (c(4)*jj-c(2)+c(5)*k2)*jj*jj+ @ (-2.d0*c(18)+(c(7)+(c(8)+c(10)*jj)*jj)*jj @ +(c(22)+(c(9)+c(11)*jj+c(12)*k2)*jj)*k2)*kl @ +(c(24)*jj+c(25)*k2+c(26)*l2)*l2*jj @ +c(27)/real(k2*k2*k2+1)*jj it=-2*iabs(mod(jmax-k-l,2))+1 if(ia.le.nad) it=-it t1(ia,na-ia+1)=it if((kk.le.j).and.(ll.le.v)) then c c...calcul de la colonne ja du premier terme non diagonal c if(ja.eq.0)then do m=-v,ll,2 if(m.eq.ll)then ib=kk else ib=jmax end if ! m do ip=-jmax,ib ir=mod(ip-m,3) if(ir.lt.0)ir=ir+3 if(ir.eq.tau)ja=ja+1 end do ! ip end do ! m end if ! ja c indla(ia)=ja ! indice colonne de q s=0.25*dsqrt(dble((v+l+2)*(v-l)*(jj-k*(k+1))* + (jj-(k+1)*kk))) deriva(ia,1)=s ! deriv q deriva(ia,2)=k2*s ! deriv qk a(ia,ja)=(c(14)+(-c(15)+c(16)*jj)*jj+c(17)*k2)*s a(ja,ia)=a(ia,ja) else indla(ia)=0 end if ! kk ja=ja+1 c case(1) ! matrice Ea ie=ie+1 ke(ie)=k le(ie)=l derive(ie,3)=-2.d0*kl ! deriv Az derive(ie,4)=l2 ! deriv X derive(ie,5)=k2 ! deriv A derive(ie,6)=k2*kl ! deriv eta.K e1(ie,ie)=c(19)*l2+(c(1)+c(13)*l2-c(3)*k2)*jj+ @ (c(20)-c(1)-c(13)*l2+(c(23)*k2-c(21)+c(6)*jj)*k2)*k2+ @ (c(4)*jj-c(2)+c(5)*k2)*jj*jj+ @ (-2.d0*c(18)+(c(7)+(c(8)+c(10)*jj)*jj)*jj @ +(c(22)+(c(9)+c(11)*jj+c(12)*k2)*jj)*k2)*kl @ +(c(24)*jj+c(25)*k2+c(26)*l2)*l2*jj @ +c(27)/real(k2*k2*k2+1)*jj c if((kk.le.j).and.(ll.le.v)) then c c calcul de la colonne ja du premier terme non diagonal c je=0 do m=-v,ll,2 if(m.eq.ll)then ib=kk else ib=jmax end if ! m do ip=-jmax,ib ir=mod(ip-m,3) if(ir.lt.0)ir=ir+3 if(ir.eq.tau)je=je+1 end do ! ip end do ! m indle(ie)=je s=0.25*dsqrt(dble((v+l+2)*(v-l)* + (jj-k*(k+1))*(jj-(k+1)*kk))) derive(ie,1)=s ! deriv q derive(ie,2)=s*k2 ! deriv qk e1(ie,je)=(c(14)+(-c(15)+c(16)*jj)*jj+c(17)*k2)*s e1(je,ie)=e1(ie,je) else indle(ie)=0 end if ! kk c case default end select ! tau c end do ! k end do ! l c if(na1.ne.na2)t1(nad+1,nad+1)=dsqrt(2.d0) do ia=1,na ! calcul de A+ et A- ˆ partir de A et T1 do ja=1,ia ai(ia,ja)=0.d0 do m=1,na s=0.d0 do n=1,na s=s+a(m,n)*t1(n,ja) end do ai(ia,ja)=ai(ia,ja)+t1(m,ia)*s end do end do end do c do ia=1,na1 do ja=1,ia a1(ia,ja)=ai(ia,ja)/2.d0 a1(ja,ia)=a1(ia,ja) end do end do c do ia=1,na2 do ja=1,ia a2(ia,ja)=ai(ia+na1,ja+na1)/2.d0 a2(ja,ia)=a2(ia,ja) end do end do c c diagonalisation, either one of the A1 or A2 matrices a2,a2 to give c eigenvalues in DA, and eigenvectors in VA c and of the E matrix to give eigenvalues in DE, eigenvectors in VE c if(na1.ge.na2) then call JACOBI(a1,na1,maxdim,da,va,nrot) ------- da(na1+1)=a2(indice,indice) else call JACOBI(a2,na2,maxdim,da,va,nrot) ------- da(na2+1)=a1(indice,indice) end if call JACOBI(e1,ne1,maxdim,de,ve,nrot) ------- c c derivatives w.r.t. q & qk c do ic=1,nq do ia=1,na do ja=1,na a(ia,ja)=0.d0 end do ! ja end do ! ia do ia=1,na ja=indla(ia) if(ja.ne.0) then a(ia,ja)=deriva(ia,ic) a(ja,ia)=deriva(ia,ic) end if end do ! ia c do ia=1,na ! A -> A1 & A2 do ja=ia,na ai(ia,ja)=0.d0 do m=1,na s=0.d0 do n=1,na s=s+a(m,n)*t1(n,ja) end do ! n ai(ia,ja)=ai(ia,ja)+t1(m,ia)*s end do ! m ai(ia,ja)=ai(ia,ja)/2.d0 ai(ja,ia)=ai(ia,ja) end do ! ja end do ! ia c do ia=1,na1 do ja=1,na1 a1(ia,ja)=ai(ia,ja) end do end do c do ia=1,na2 do ja=1,na2 a2(ia,ja)=ai(ia+na1,ja+na1) end do end do c if(na1.ge.na2) then do ia=1,na1 deriva(ia,ic)=0.d0 do m=1,na1 s=0.d0 do n=1,na1 s=s+a1(m,n)*va(n,ia) end do ! n deriva(ia,ic)=deriva(ia,ic)+va(m,ia)*s end do ! m end do ! ia deriva(na1+1,ic)=a2(indice,indice) else ! na1|na2 do ia=1,na2 deriva(ia,ic)=0.d0 do m=1,na2 s=0.d0 do n=1,na2 s=s+a2(m,n)*va(n,ia) end do ! n deriva(ia,ic)=deriva(ia,ic)+va(m,ia)*s end do ! m end do ! ia deriva(na2+1,ic)=a1(indice,indice) end if ! na1|na2 c do ie=1,ne1 do je=1,ne1 e1(ie,je)=0.d0 end do ! je end do ! ie c do ie=1,ne1 je=indle(ie) if(je.ne.0) then e1(ie,je)=derive(ie,ic) e1(je,ie)=derive(ie,ic) end if end do ! ie c do ie=1,ne1 derive(ie,ic)=0.d0 do m=1,ne1 s=0.d0 do n=1,ne1 s=s+e1(m,n)*ve(n,ie) end do ! n derive(ie,ic)=derive(ie,ic)+ve(m,ie)*s end do ! m end do ! ie end do ! ic c c derivatives w.r.t. Az, X, A c do ic=3,nx do ie=1,ne1 do je=1,ne1 e1(ie,je)=ve(ie,je)*derive(ie,ic) end do ! je end do ! ie do ie=1,ne1 derive(ie,ic)=0.d0 do je=1,ne1 derive(ie,ic)=derive(ie,ic)+ve(je,ie)*e1(je,ie) end do ! je end do ! ie c do ia=1,na do ja=1,ia ai(ia,ja)=0.d0 do m=1,na ai(ia,ja)=ai(ia,ja)+t1(m,ia)*t1(m,ja)*deriva(m,ic) end do ! m ai(ia,ja)=ai(ia,ja)/2.d0 ai(ja,ia)=ai(ia,ja) end do ! ja end do ! ia do ia=1,na1 do ja=1,na1 a1(ia,ja)=ai(ia,ja) end do end do c do ia=1,na2 do ja=1,na2 a2(ia,ja)=ai(ia+na1,ja+na1) end do end do c if(na1.ge.na2) then do ia=1,na1 deriva(ia,ic)=0.d0 do m=1,na1 s=0.d0 do n=1,na1 s=s+a1(m,n)*va(n,ia) end do ! n deriva(ia,ic)=deriva(ia,ic)+va(m,ia)*s end do ! m end do ! ia deriva(na1+1,ic)=a2(indice,indice) else ! na1|na2 do ia=1,na2 deriva(ia,ic)=0.d0 do m=1,na2 s=0.d0 do n=1,na2 s=s+a2(m,n)*va(n,ia) end do ! n deriva(ia,ic)=deriva(ia,ic)+va(m,ia)*s end do ! m end do ! ia deriva(na2+1,ic)=a1(indice,indice) end if ! na1|na2 end do ! ic c c stockage c do ic=i,numl jm=2*(j-jn(ic))-1 if((jm.ne.-1).and.(jm.ne.1)) cycle kle=kn(ic)*ln(ic) if(mod(kn(ic)-ln(ic),3).eq.0) then do ia=1,max0(na1+1,na2+1) if(kle.ne.ka(ia)*la(ia)) cycle if(ln(ic).eq.0) then ! PAS BON POUR V=3 if(iabs(kn(ic)).ne.iabs(ka(ia))) cycle else if(((kle.eq.1).and.(v.eq.1)).or.((kle.eq.4).and.(v.eq.2))) then ir=1 if(jm.eq.1)ir=-1 ! doublet k if(kn(ic).ne.ka(ia)*ir) cycle if(ln(ic).ne.la(ia)*ir) cycle end if ! kle=1;4 end if ! la=0 fber(ic)=fber(ic)+jm*da(ia) amat(ic,14)=amat(ic,14)+jm*deriva(ia,1) amat(ic,15)=amat(ic,15)-jm*jj*deriva(ia,1) amat(ic,16)=amat(ic,16)+jm*jj*jj*deriva(ia,1) amat(ic,17)=amat(ic,17)+jm*deriva(ia,2) amat(ic,18)=amat(ic,18)+jm*deriva(ia,3) amat(ic,19)=amat(ic,19)+jm*deriva(ia,4) amat(ic,20)=amat(ic,20)+jm*deriva(ia,5) amat(ic,22)=amat(ic,22)+jm*deriva(ia,6) exit end do ! ia else do ie=1,ne1 if(kle.ne.ke(ie)*le(ie)) cycle if(kle.eq.0) then if(iabs(kn(ic)).ne.iabs(ke(ie))) cycle if(iabs(ln(ic)).ne.iabs(le(ie))) cycle end if ! le=0 fber(ic)=fber(ic)+jm*de(ie) amat(ic,14)=amat(ic,14)+jm*derive(ie,1) amat(ic,15)=amat(ic,15)-jm*jj*derive(ie,1) amat(ic,16)=amat(ic,16)+jm*jj*jj*derive(ie,1) amat(ic,17)=amat(ic,17)+jm*derive(ie,2) amat(ic,18)=amat(ic,18)+jm*derive(ie,3) amat(ic,19)=amat(ic,19)+jm*derive(ie,4) amat(ic,20)=amat(ic,20)+jm*derive(ie,5) amat(ic,22)=amat(ic,22)+jm*derive(ie,6) exit end do ! ie end if ! k-l=3n end do ! ic end if ! lg end do ! in end do ! i c c...Calculate transition frequencies c if(lg4)then write(1,1010) else write(1,2002) endif c do i=1,numl j=jn(i) c c...derivatives derivatives which can be evaluated analytically - note c that these are now derivatives in the transition frequency c if(j.lt.0) then ! doublet l, v=1 j=-j jj=j*(j+1) amat(i,14)=jj amat(i,15)=-jj*jj amat(i,16)=-amat(i,15)*jj amat(i,17)=jj fber(i)=amat(i,14)*c(14)+amat(i,15)*c(15)+amat(i,16)*c(16) @ +amat(i,17)*c(17) else jj=j+1 k2=kn(i)*kn(i) kl=kn(i)*ln(i) amat(i,1)=2*jj amat(i,2)=-4*jj**3 amat(i,3)=-2*jj*k2 amat(i,4)=2*jj**3*dble(3*j*j+6*j+4) amat(i,5)=-amat(i,2)*k2 amat(i,6)=-amat(i,3)*k2 amat(i,7)=2*jj*kl amat(i,8)=-amat(i,2)*kl amat(i,9)=-amat(i,3)*kl amat(i,10)=2*jj**3*dble(3*j*j+6*j+4)*kl amat(i,11)=amat(i,8)*k2 amat(i,12)=amat(i,9)*k2 amat(i,13)=2*jj*ln(i)*ln(i) amat(i,24)=4*jj**3*ln(i)*ln(i) amat(i,25)=2*jj*k2*ln(i)*ln(i) amat(i,26)=2*jj*ln(i)*ln(i)*ln(i)*ln(i) amat(i,27)=real(2*jj)/real(k2*k2*k2+1) end if ! j<0 c verb(i)=fber(i)-fr(i) h1=w(i) h2=verb(i)*h1 c kmls=kn(i)-ln(i) if(ln(i).lt.0)kmls=-kmls c if(lg4) then hat=sdev(i) ec=-h2/dsqrt(varz*(1-hat)) ! ti if(abs(ec).gt.3.5) then testi='***' else testi=' ' end if ! ec>3.5 pp=f2/w(i)*1000.d0 write(1,1011)i,jn(i),kn(i),ln(i),fr(i),-verb(i),ec, + testi,hat,pp,fber(i),kmls else ! lg4 if((kn(i)/3)*3.eq.kn(i))then gi=(1.0+2.0/(2.0*spnucl+1.0)**2) else gi=(1.0-2.0/(2.0*spnucl+1.0)**2) endif if(kn(i).eq.0)then gk=1.0 else gk=2.0 endif rint(i)=gi*gk*(1.0-real(kn(i)**2)/real((jn(i)+1)**2)) if(kn(i).eq.ln(i))rint(i)=rint(i)*0.5 pp=f2/w(i)*1000.d0 write(1,1009)i,jn(i),kn(i),ln(i),fr(i),-verb(i),fber(i), + kmls,rint(i),pp end if ! lg4 c sdev(i)=h2 dacc=sdev(i)**2+dacc do ib=1,nconst amat(i,ib)=amat(i,ib)*h1 end do ! ib end do ! i c c...sort accoring to frequency c if(.not.lg4)then do 2004 isort=1,numl wk(isort)=fber(isort) ipt(isort)=isort 2004 continue call sorth(1,numl) ------- write(1,'(1x/'' Lines sorted in order of calculated'', * '' frequency:''/)') do 2005 isort=1,numl i=ipt(isort) kmls=kn(i)-ln(i) if(ln(i).lt.0)kmls=-kmls pp=f2/w(i)*1000.d0 write(1,1009)i,jn(i),kn(i),ln(i),fr(i),-verb(i), * fber(i), kmls,rint(i),pp 2005 continue endif c c...sum of squared defects c c write(6,1012)dacc write(1,1012)dacc c c...MM is the number of fitted constants, compress AMAT to eliminate fixed c vectors c mm=nconst if(.not.lg4) then do ic=1,icst mm=mm-1 n=icc(ic)-ic+1 if(n.gt.mm)cycle do i=n,mm do l=1,numl amat(l,i)=amat(l,i+1) end do ! l end do ! i end do ! ic end if ! lg4 c c...standard deviation now printed in units of frequency (zk) c dacc=dsqrt(dacc/dble(numl-mm))*1000.d0 write(6,1013)dacc/1000.d0 write(1,1013)dacc/1000.d0 c c...write new header for data file on exit c if(lg4) then open(3,file='v2e.hdr',status='unknown') do ic=1,nconst c write(6,1006)c(ic),ilib(ic),param(ic) write(3,1006)c(ic),ilib(ic),param(ic) end do ! ic close(3) write(1,'(1x,78(1h-))') close(1) stop ! exit from program end if c write(6,1016)itm write(1,1016)itm c c...DA = -change to constant value, Y =standard deviation in constant c mm= number of fitted constants c call FITSVD(sdev,da,y,maxlin,maxdim,numl,mm,.true.) ------- varz=amat(1,mm+1) c c... expand DA so that it corresponds to C, ie. zero values are inserted c in DA(i) if i'th constant is fixed c do ic=1,icst n=icc(ic) if(n.le.mm)then do l=n,mm i=mm+n-l da(i+1)=da(i) end do ! l end if ! n c if(n.gt.22)exit ! HK pas fittable da(n)=0.d0 mm=mm+1 end do ! ic c c...evaluation of corrected constants c do i=1,nconst c(i)=c(i)-da(i) end do ! i c end do ! itm (end of main fitting loop) c c . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . c 1001 format(a6,a1,f25.13,2f22.13) 1005 format(6h NUML ,i3,3x,4hNORM,d10.3,3x,4hHAT0,e10.3/) 1006 format(g18.10,i4,3x,a6) 1009 format(1h ,i3,'*',3i3,f12.3,f10.3,f12.3,i5,f12.4,f12.2) 1010 format(1x/7x,'J K l Fobs o-c t(i) hii ', * ' error Fcalc'/) 2002 format(1x/7x,'J K l Fobs o-c Fcalc ', * '(K-l)sgn(l) I.rel error'/) 1011 format(1h ,i3,'*',3i3,f12.3,f10.3,f7.2,a3,f7.2,f10.1,f12.3,i5) 1012 format(/23h sum of squared defects,f15.6) 1013 format(/23h standard deviation,f15.6) 1016 format(1x/' =============> ITERATION number',i3,1x,40(1h-)/) c stop end c c****************************************************************************** C SUBROUTINE JACOBI(A,N,NP,D,V,NROT) C ***Numerical Recipes p. 346*** c c Eigenvalues and eigenvectors of A, which is partly destroyed in c the process c A - array A(N,N) to be diagonalized stored in physical array c A(NP,NP) c D - eigenvalues on output c V - eigenvectors on output c NROT - number of rotations required returned on output c c IMPLICIT DOUBLE PRECISION (A-H,P-Z) PARAMETER (MAXDIM=100) DIMENSION A(NP,NP),D(NP),V(NP,NP),B(MAXdim),Z(MAXdim) c DO IP=1,N DO IQ=1,N V(IP,IQ)=0. end do V(IP,IP)=1. end do DO IP=1,N B(IP)=A(IP,IP) D(IP)=B(IP) Z(IP)=0. end do NROT=0 c DO I=1,50 SM=0. DO IP=1,N-1 DO IQ=IP+1,N SM=SM+DABS(A(IP,IQ)) end do end do IF(SM.EQ.0.) RETURN c c...here correction relative to version from J.DEMAISON: from IF(L.LT.. c to IF(I.LT.....) consistent with listing in NUM.REC. (zk) c IF(I.LT.4) THEN TRESH=0.2*SM/N**2 ELSE TRESH=0. END IF DO IP=1,N-1 DO IQ=IP+1,N G=100.*DABS(A(IP,IQ)) IF((I.GT.4).AND.(DABS(D(IP))+G.EQ.DABS(D(IP))) * .AND.(DABS(D(IQ))+G.EQ.DABS(D(IQ)))) THEN A(IP,IQ)=0. ELSE IF(DABS(A(IP,IQ)).GT.TRESH) THEN H=D(IQ)-D(IP) IF(DABS(H)+G.EQ.DABS(H)) THEN T=A(IP,IQ)/H ELSE THETA=0.5*H/A(IP,IQ) T=1./(DABS(THETA)+DSQRT(1.+THETA**2)) IF(THETA.LT.0.)T=-T END IF C=1./DSQRT(1+T**2) S=T*C TAU=S/(1.+C) H=T*A(IP,IQ) Z(IP)=Z(IP)-H Z(IQ)=Z(IQ)+H D(IP)=D(IP)-H D(IQ)=D(IQ)+H A(IP,IQ)=0. DO J=1,IP-1 G=A(J,IP) H=A(J,IQ) A(J,IP)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) end do DO J=IP+1,IQ-1 G=A(IP,J) H=A(J,IQ) A(IP,J)=G-S*(H+G*TAU) A(J,IQ)=H+S*(G-H*TAU) end do DO J=IQ+1,N G=A(IP,J) H=A(IQ,J) A(IP,J)=G-S*(H+G*TAU) A(IQ,J)=H+S*(G-H*TAU) end do DO J=1,N G=V(J,IP) H=V(J,IQ) V(J,IP)=G-S*(H+G*TAU) V(J,IQ)=H+S*(G-H*TAU) end do NROT=NROT+1 END IF end do end do DO IP=1,N B(IP)=B(IP)+Z(IP) D(IP)=B(IP) Z(IP)=0. end do end do c PAUSE '50 iterations should never happen' c RETURN END C C**************************** SVD ******************************************** C SUBROUTINE SVD(M,N,Q) c IMPLICIT INTEGER*4 (I-N) parameter (nconst=27,maxlin=330,maxdim=100) REAL*8 U,V,Q,E,C,F,G,H,S,X,Y,Z,EPS,TOL,cdummy DIMENSION U(maxlin,nconst),V(maxdim,maxdim),Q(maxdim), + E(maxlin),cdummy(maxdim,maxdim) common /work/v,cdummy COMMON /DERIV/U c LOGICAL LGU,LGV,LG LGU=.FALSE. LGV=.FALSE. EPS=1.0D-15 TOL=2.3D-308 ! 0.22D-2920 extended G=0.D0 X=0.D0 DO 2 I=1,N E(I)=G S=0.D0 L=I+1 LG=L.EQ.N+1 DO 3 J=I,M S=U(J,I)**2+S 3 continue G=0.D0 IF(S-TOL)9,5,5 5 F=U(I,I) c c...in original program extended precision routine QSQRT was used (zk) G=DSQRT(S) IF(F.GE.0.)G=-G H=F*G-S U(I,I)=F-G IF(LG)GOTO 9 DO 6 J=L,N S=0.D0 DO 7 K=I,M S=U(K,I)*U(K,J)+S 7 continue F=S/H DO 8 K=I,M U(K,J)=U(K,I)*F+U(K,J) 8 continue 6 CONTINUE 9 Q(I)=G IF(LG)GOTO 4 S=0.D0 DO 10 J=L,N S=U(I,J)**2+S 10 continue G=0.D0 IF(S-TOL)4,11,11 11 F=U(I,I+1) G=DSQRT(S) IF(F.GE.0.)G=-G H=F*G-S U(I,I+1)=F-G DO 12 J=L,N 12 E(J)=U(I,J)/H DO 13 J=L,M S=0. DO 14 K=L,N 14 S=U(J,K)*U(I,K)+S DO 13 K=L,N 13 U(J,K)=E(K)*S+U(J,K) 4 Y=DABS(Q(I))+DABS(E(I)) IF(Y.GT.X)X=Y 2 CONTINUE IF(LGV)GOTO 15 DO 16 II=1,N I=N+1-II LG=L.EQ.N+1 IF(LG)GOTO 211 IF(G.EQ.0.)GOTO 17 H=U(I,I+1)*G DO 18 J=L,N 18 V(J,I)=U(I,J)/H DO 19 J=L,N S=0.D0 DO 20 K=L,N 20 S=U(I,K)*V(K,J)+S DO 19 K=L,N 19 V(K,J)=V(K,I)*S+V(K,J) 17 DO 21 J=L,N V(I,J)=0.D0 21 V(J,I)=0.D0 211 V(I,I)=1.D0 G=E(I) 16 L=I 15 IF(LGU)GOTO 22 DO 23 II=1,N I=N+1-II L=I+1 G=Q(I) LG=L.EQ.N+1 IF(LG)GOTO 241 DO 24 J=L,N U(I,J)=0.D0 24 continue 241 IF(G.EQ.0.)GOTO 25 H=U(I,I)*G IF(LG)GOTO 261 DO 26 J=L,N S=0.D0 DO 27 K=L,M S=U(K,I)*U(K,J)+S 27 continue F=S/H DO 26 K=I,M U(K,J)=U(K,I)*F+U(K,J) 26 continue 261 DO 28 J=I,M U(J,I)=U(J,I)/G 28 continue GOTO 23 25 DO 29 J=I,M U(J,I)=0.D0 29 continue 23 U(I,I)=U(I,I)+1. 22 EPS=EPS*X c DO 30 II=1,N K=N+1-II 31 DO 32 LL=1,K L=K+1-LL IF(DABS(E(L)).LE.EPS)GOTO 33 IF(DABS(Q(L-1)).LE.EPS)GOTO 34 32 CONTINUE 34 C=0.D0 S=1.D0 L1=L-1 DO 35 I=L,K F=E(I)*S E(I)=E(I)*C IF(DABS(F).LE.EPS)GOTO 33 G=Q(I) H=DSQRT(F*F+G*G) Q(I)=H C=G/H S=-F/H IF(LGU)GOTO 35 DO 36 J=1,M Y=U(J,L1) Z=U(J,I) U(J,L1)=Y*C+Z*S U(J,I)=-Y*S+Z*C 36 continue 35 CONTINUE 33 Z=Q(K) IF(L.EQ.K)GOTO 37 X=Q(L) Y=Q(K-1) G=E(K-1) H=E(K) F=H*Y*2. F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/F G=DSQRT(F*F+1.) S=F-G IF(F.GE.0.)S=F+G F=((X-Z)*(X+Z)+(Y/S-H)*H)/X C=1.D0 S=1.D0 LL=L+1 DO 38 I=LL,K G=E(I) Y=Q(I) H=S*G G=C*G Z=DSQRT(F*F+H*H) E(I-1)=Z C=F/Z S=H/Z F=X*C+G*S G=-X*S+G*C H=Y*S Y=Y*C IF(LGV)GOTO 39 DO 40 J=1,N X=V(J,I-1) Z=V(J,I) V(J,I-1)=X*C+Z*S V(J,I)=-X*S+Z*C 40 continue 39 Z=DSQRT(F*F+H*H) Q(I-1)=Z C=F/Z S=H/Z F=C*G+S*Y X=-S*G+C*Y IF(LGU)GOTO 38 DO 41 J=1,M Y=U(J,I-1) Z=U(J,I) U(J,I-1)=Y*C+Z*S U(J,I)=-Y*S+Z*C 41 continue 38 CONTINUE E(L)=0.D0 E(K)=F Q(K)=X GOTO 31 37 IF(Z.GE.0.)GOTO 30 Q(K)=-Z IF(LGV)GOTO 30 DO 42 J=1,N V(J,K)=-V(J,K) 42 continue 30 CONTINUE c RETURN END C C******************************************* FITSVD *********************** C SUBROUTINE FITSVD(B,X,Y,MA,NA,M,N,LG) c c Related to routine SVDFIT from NumRec p.518 c c MA - total number of possible data points c NA - total number of coefficients c c M - number of data points c N - number of coeefficients to be determined c X - corrections to constants on output c Y - standard deviations on constants on output c c C - correlation matrix on output c V - eigenvectors c IMPLICIT INTEGER*4 (I-N) PARAMETER (nconst=27, maxlin=330, maxdim=100) REAL*8 A,B,X,V,F,Y,Z,C,BN,S,T,U DIMENSION A(maxlin,nconst),V(maxdim,maxdim),B(maxlin),X(maxdim), + C(maxdim,maxdim),F(nconst),Z(nconst),Y(NA) REAL*4 PHI(nconst),PI(nconst,nconst),RES(maxlin),W CHARACTER*1 CHDF(nconst),CHTI,CHH,CHD COMMON /DERIV/A common /work/v,c LOGICAL LG,LGU c 1000 FORMAT(40H NORMALIZATION OF EQUATIONS OF CONDITION) 1001 FORMAT(10X,D15.7,E14.3) 1002 FORMAT(I3,20F7.3) 1003 FORMAT(20D13.5) ! format de eigenvalues 1004 FORMAT(17H CONDITION NUMBER,D14.4) C 1005 FORMAT(12H DETERMINANT,5X,D14.4) 1006 FORMAT(21H NEAR SINGULAR MATRIX,D14.4) C 1007 FORMAT(5H NORM,12X,D14.4) C 1008 FORMAT(9H DET/NORM,8X,D14.4) 1009 FORMAT(1x/9H VARIANCE,8X,f15.6) !/18H SUM DEFECT VECTOR) C 1010 FORMAT(17H REDHEFFER ERROR ,E14.4) 1011 FORMAT(' CHANGE IN CONSTANT, STANDARD ERROR') 1012 FORMAT(2X,2F25.13) 1013 FORMAT(1x/19H CORRELATION MATRIX) 1014 FORMAT(19F7.3) 1015 FORMAT(27H VARIANCE-COVARIANCE MATRIX) 1016 FORMAT(5X,9D14.6) 1017 FORMAT(I8,12F6.3,(8X,12F6.3)) 1042 FORMAT(1x/' Leverage Components uij^2, threshold:',F8.3) c LGU=.FALSE. C LG=.TRUE. detailed output and tests C LGU=.TRUE. Leverage Components uij^2 L=6 ! n¡ output file c c IF(LG) WRITE(L,1000) IF(LG) WRITE(1,1000) DO 1 I=1,N S=0.d0 DO 2 J=1,M S=A(J,I)**2+S 2 CONTINUE if(s.eq.0.d0)then s=1.d100 else S=1.D0/DSQRT(S) endif c IF(LG) WRITE(L,1001) S IF(LG) WRITE(1,1001) S F(I)=S DO 1 J=1,M A(J,I)=A(J,I)*S 1 CONTINUE S=0.d0 DO 3 I=1,M RES(I)=B(I) S=B(I)**2+S 3 CONTINUE BN=1.D0/DSQRT(S) c IF(LG) WRITE(L,1001) BN IF(LG) WRITE(1,1001) BN DO 4 I=1,M 4 B(I)=B(I)*BN CALL SVD(M,N,X) ------- IF(LG) THEN c WRITE(L,*)'SINGULAR VALUES' WRITE(1,'(1x/1x,a)')'SINGULAR VALUES' c WRITE(L,1003) (X(I),I=1,N) WRITE(1,1003) (X(I),I=1,N) c WRITE(L,*)'Eigenvectors V' WRITE(1,'(1x/1x,a)')'Eigenvectors V' DO I=1,N c WRITE(L,1002)I,(V(I,J),J=1,N) ! eigenvectors WRITE(1,1002)I,(V(I,J),J=1,N) ! eigenvectors END DO THRES=FLOAT(N)/FLOAT(M) IF(LGU) THEN c WRITE(L,1042)THRES WRITE(1,1042)THRES DO I=1,M c WRITE(L,1002)I,(A(I,J)**2,J=1,N) ! UIJ^2 WRITE(1,1002)I,(A(I,J)**2,J=1,N) ! UIJ^2 END DO END IF END IF S=X(1) T=S DO 6 I=2,N S=DMAX1(S,X(I)) SMAX=S T=DMIN1(T,X(I)) 6 continue IF(T) 7,8,7 8 WRITE(L,1006) T WRITE(1,1006) T STOP 7 T=S/T IF(T.GT.1D16) GOTO 8 c IF(LG) WRITE(L,1004) T IF(LG) WRITE(1,1004) T DO 9 I=1,N T=0 DO 10 J=1,M T=A(J,I)*B(J)+T 10 continue Y(I)=T S=X(I) DO 9 J=1,N C(I,J)=V(J,I)*S 9 continue T=0.D0 DO 11 I=1,N S=0.D0 DO 12 J=1,N V(I,J)=V(I,J)/X(J) S=C(J,I)*Y(J)+S 12 continue Z(I)=S 11 CONTINUE DO I=1,N PHI(I)=0. DO J=1,N PI(I,J)=V(I,J)*V(I,J) PHI(I)=PHI(I)+PI(I,J) END DO DO J=1,N PI(I,J)=PI(I,J)/PHI(I) END DO END DO c IF(LG) THEN c WRITE(L,*)'Variance-Decomposition Proportions' WRITE(1,'(1x/1x,a)')'Variance-Decomposition Proportions' c WRITE(L,*)'Cond in v(b1) v(b2) v(b3) v(b4) v(b5) v(b6)' WRITE(1,*)'Cond in v(b1) v(b2) v(b3) v(b4) v(b5) v(b6)' DO I=1,N c WRITE(L,1017)NINT(SMAX/X(I)),(PI(J,I),J=1,N) WRITE(1,1017)NINT(SMAX/X(I)),(PI(J,I),J=1,N) END DO END IF DO 14 I=1,N S=0.D0 DO 15 J=1,N S=V(I,J)*Y(J)+S 15 continue X(I)=S 14 continue DO 27 I=1,M DO 28 J=1,N Y(J)=A(I,J) 28 continue DO 27 J=1,N U=1.D0/F(J) S=0.D0 DO 29 K=1,N S=C(K,J)*Y(K)+S 29 continue A(I,J)=S*U 27 continue U=0.D0 DO 16 I=1,N U=X(I)*Z(I)+U S=0.D0 DO 17 J=1,N S=C(I,J)*X(J)+S 17 continue Y(I)=S 16 continue U=(1.D0-U)/(DFLOAT(M-N)*BN*BN) c IF(LG) WRITE(L,1009) U IF(LG) WRITE(1,1009) U A(1,N+1)=U ! variance VARZ=U DO 20 I=1,N DO 20 J=1,I S=0.D0 DO 21 K=1,N S=V(I,K)*V(J,K)+S 21 CONTINUE C(I,J)=S C(J,I)=S 20 CONTINUE T=DSQRT(U) c IF(LG) WRITE(L,1011) c IF(LG) WRITE(1,1011) DO 22 I=1,N X(I)=X(I)*F(I)/BN ! solution ie. B.i-B.i+1 or -change Z(I)=DSQRT(C(I,I)) Y(I)=Z(I)*F(I)*T ! standard deviation of parameter I c IF(LG) WRITE(L,1012) -X(I),Y(I) c IF(LG) WRITE(1,1012) -X(I),Y(I) 22 CONTINUE c c...correlation matrix c c IF(LG) WRITE(L,1013) IF(LG) WRITE(1,1013) DO 23 I=1,N S=Z(I) DO 24 J=1,N T=Z(J)*S C(I,J)=C(I,J)/T 24 CONTINUE c IF(LG) WRITE(L,1014)(C(I,J),J=1,I) IF(LG) WRITE(1,1014)(C(I,J),J=1,I) 23 CONTINUE c IF(LGU) WRITE(L,1015) ! matrice des variance-covariances IF(LGU) WRITE(1,1015) ! matrice des variance-covariances DO 25 I=1,N S=Y(I) DO 26 J=1,N V(I,J)=C(I,J)*Y(J)*S 26 CONTINUE c IF(LGU) WRITE(L,1016)(V(I,J),J=1,I) ! matrice des variance-covariances IF(LGU) WRITE(1,1016)(V(I,J),J=1,I) ! matrice des variance-covariances 25 CONTINUE THRESDF=2/SQRT(FLOAT(M)) THRESH=2*FLOAT(N)/FLOAT(M) THRESD=2*SQRT(FLOAT(N)/FLOAT(M-N)) c c...statistical tests c C IF(LG)WRITE(L,1041)THRESH,THRESD,THRESDF IF(LG)WRITE(1,1041)THRESH,THRESD,THRESDF 1041 FORMAT(1h /' TESTS, thresholds: hii',F6.3,' DFFITS',F6.3 @,' CD2 1 DFBETAs',F6.3/3x,'i',2x,'rsd',3x,'t(i)',3x,'hii' @,3x,'DFFITS',2X,'CD2',5X,'DFBETAS') c DO 30 I=1,M S=0.D0 W=0 DO J=1,N U=A(I,J) c W=W+FLOAT(U*X(J)) original W=W+ U*X(J) T=0 DO 31 K=1,N BN=V(K,J)*A(I,K)/VARZ T=T+BN S=BN*U+S 31 continue F(J)=T END DO B(I)=S ! hii IF(LG) THEN CHH=" " CHTI=" " CHD=" " RES(I)=RES(I)-W ! residual IF(S.GT.THRESH) CHH="*" VARZI=(FLOAT(M-N)*VARZ-RES(I)**2/(1-S))/FLOAT(M-N-1) ! s(i)^2 TI=ABS(RES(I)/SQRT((1-S)*VARZI)) ! t(i) IF(TI.GT.3.5) CHTI="*" DFFIT=RES(I)*SQRT(S/VARZI)/(1-S) IF(ABS(DFFIT).GT.THRESD) CHD="*" CD2=RES(I)**2*S/(FLOAT(N)*VARZ*(1-S)**2) ! Cook 1043 FORMAT(I4,F6.2,3(F6.2,A1),F6.2,2x,15(F6.2,A1)) DO J=1,N CHDF(J)=" " F(J)=F(J)*RES(I)/((1-B(I))*Y(J)) ! DFBETAS IF(ABS(F(J)).GT.THRESDF) CHDF(J)="*" END DO C WRITE(L,1043)I,RES(I)/SQRT(VARZ),TI,CHTI,S,CHH,DFFIT,CHD,CD2 C @,(F(J),CHDF(J),J=1,N) WRITE(1,1043)I,RES(I)/SQRT(VARZ),TI,CHTI,S,CHH,DFFIT,CHD,CD2 @,(F(J),CHDF(J),J=1,N) END IF 30 CONTINUE c RETURN END C C_____________________________________________________________________________ C SUBROUTINE SORTH(NSTART,N) c c This routine is based on the SORT2 'heapsort' routine from Numerical c Recipes and sorts the quantities in vector WK from WK(NSTART) to WK(N) C in ascending order of magnitude and also accordingly rearranges vector C IPT of pointers to original positions of sorted quantities. c PARAMETER (MAXLIN=330) c COMMON /sort/WK,ipt INTEGER*2 IPT(MAXLIN),IIPT,L,I,J,IR REAL*8 WK(MAXLIN),WWK C L=N/2+1 IR=N 10 CONTINUE IF(L.GT.NSTART)THEN L=L-1 WWK=WK(L) IIPT=IPT(L) ELSE WWK=WK(IR) IIPT=IPT(IR) WK(IR)=WK(1) IPT(IR)=IPT(1) IR=IR-1 IF(IR.EQ.NSTART)THEN WK(1)=WWK IPT(1)=IIPT RETURN ENDIF ENDIF I=L J=L+L 20 IF(J.LE.IR)THEN IF(J.LT.IR)THEN IF(WK(J).LT.WK(J+1))J=J+1 ENDIF IF(WWK.LT.WK(J))THEN WK(I)=WK(J) IPT(I)=IPT(J) I=J J=J+J ELSE J=IR+1 ENDIF GO TO 20 ENDIF WK(I)=WWK IPT(I)=IIPT GO TO 10 c RETURN END C C_____________________________________________________________________________ C_____________________________________________________________________________