; TO DO :
;
;
;
; $Id: g98log__define.pro,v 1.3 2005/08/19 21:49:37 craig_nist Exp $
;+
; NAME:
;       G98LOG__DEFINE
;
; PURPOSE:
;
;       Object class for Gaussian98/03 output file
;
; AUTHORS:
;       Craig Brown
;     NIST Center for Neutron Research
;       100 Bureau Drive
;     Gaithersburg, MD 20899
;       Phone: (301) 975-5134
;       E-mail: craig.brown@nist.gov
;
; CATEGORY:
;
;       Objects, widgets, data interpretation
;
; CALLING SEQUENCE:
;
;       object = obj_new('G98LOG',DATADIRECTORY = datadirectory)
;
;
;###############################################################################
; LICENSE:
;  The software in this file is written by an employee of
;  National Institute of Standards and Technology
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################
;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function gausseditatoms,group_leader=base,in,out,changed
gaussvarEdit,out,in,changed,group=base,y_scroll_size=(n_elements(in.list))
return,out
return,changed
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98log::cleanup
;remove pointers
if ptr_valid(self.isotopeinfoptr) then ptr_free,self.isotopeinfoptr
if ptr_valid(self.initialorientationptr) then ptr_free,self.initialorientationptr
if ptr_valid(self.atomorigptr) then ptr_free,self.atomorigptr
if ptr_valid(self.atomnewptr) then ptr_free, self.atomnewptr
if ptr_valid(self.symmodeptr) then ptr_free,self.symmodeptr
if ptr_valid(self.origfreqptr) then ptr_free,self.origfreqptr
if ptr_valid(self.freqptr) then ptr_free,self.freqptr
if ptr_valid(self.redmassptr) then ptr_free,self.redmassptr
if ptr_valid(self.frccptr) then ptr_free,self.frccptr
if ptr_valid(self.iriptr) then ptr_free,self.iriptr
if ptr_valid(self.ramanptr) then ptr_free,self.ramanptr
if ptr_valid(self.depolarptr) then ptr_free,self.depolarptr
if ptr_valid(self.vectorptr) then ptr_free,self.vectorptr
if ptr_valid(self.eigenyptr) then ptr_free,self.eigenyptr
if ptr_valid(self.massevecsptr) then ptr_free,self.massevecsptr
if ptr_valid(self.spectrumptr) then ptr_free,self.spectrumptr


return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98log::quit
widget_control,self.tlb,/destroy
return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function g98log::edit_atoms,Group_Leader=Group_leader

              list=(*self.atomorigptr).atomlist
              atom=(*self.atomnewptr).atomlist
              ixsection=(*self.atomorigptr).ixsection
              IXS=(*self.atomnewptr).ixsection
              mass=(*self.atomnewptr).mass
              origmass=(*self.atomorigptr).mass

              inatom={list:list,$
              origmass:origmass,$
              ixsection:ixsection}

              outatom={atom:atom,$
              mass:mass,$
              IXS:IXS}
              changed=0

              atli=gausseditatoms(group_leader=Group_leader,inatom,outatom,changed)

              if (changed)  then begin
              (*self.atomnewptr).atomlist=atli.atom
              (*self.atomnewptr).mass=atli.mass
              (*self.atomnewptr).ixsection=atli.IXS
              ret=self->eigenvec()
              return,1
              endif
return,0
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

function g98log::reduce,mono=mono,coll=coll,scalefactor=scalefactor,$
              oversf=oversf,Freqfactor=Freqfactor,$
              flatbg=flatbg,linearbg=linearbg,DWF=DWF,plotover=plotover
compile_opt hidden
;
;Purpose: On some platforms, when IDLgrWindow::Draw is invoked, math
;errors (e.g. "% Program caused arithmetic error: Floating illegal
;operand") are printed.  xobjview::draw suppresses the
;printing of these errors.
;
;Flush and print any accumulated math errors
;
void = check_math(/print)
;
;Silently accumulate any subsequent math errors.
;
orig_except = !except
!except = 0


if keyword_set(linearbg) then b1=linearbg else b1=0
if keyword_set(flatbg) then b0=flatbg else b0=0
if keyword_set(plotover) then plotover=plotover else plotover=0

treatment=['Filename: ' +self.filename]
treatment=[[treatment],'Gaussian job run: '+strjoin(SELF.job)]
treatment=[[treatment],'max. force= '+string((SELF.maxforce)[0])]
treatment=[[treatment],'rms. force= '+string((SELF.rmsforce)[0])]
treatment=[[treatment],'max. disp= '+string((SELF.maxdisp)[0])]
treatment=[[treatment],'rms. disp= '+string((SELF.rmsdisp)[0])]
treatment=[[treatment],'numberatoms= '+string(SELF.numberatoms)]
;____________________________________________________________________________________
;START RRESOLTION PART
d_cu=1.278                 ; Cu(220) d-spacing A
d_pg=3.354
dmono=d_cu
coll1=20.0
coll2=20.0
fundamental=fltarr(2901)
overtone=fltarr(2901)
plotspec=fltarr(2901)
freqspec=findgen(2901)
freqspec=((freqspec/10)+10)*8.06544

case mono of
    0: begin
           dmono=d_cu
       end
    1: begin
       dmono=d_pg
    end
    2: begin
       ;gaussplotcalc,event
       return,0
    end
endcase
case coll of
    0: begin
           coll1=10.0
           coll2=10.0
       end
    1: begin
        coll1=10.0
        coll2=20.0
       end
    2: begin
        coll1=20.0
        coll2=10.0
       end
    3: begin
        coll1=20.0
        coll2=20.0
       end
    4: begin
        coll1=20.0
        coll2=40.0
       end
    5: begin
        coll1=40.0
        coll2=20.0
       end
    6: begin
        coll1=40.0
        coll2=40.0
       end
    7: begin
        coll1=40.0
        coll2=60.0
       end
    8: begin
        coll1=60.0
        coll2=40.0
       end
    9: begin
        coll1=60.0
        coll2=60.0
       end
endcase

col1=coll1/60.0/57.269
col2=coll2/60.0/57.269
rmos=30./60./57.269
tmpd=col1^2+col2^2+4.*rmos^2
tmpn=(col1*col2)^2+(col1*rmos)^2+(col2*rmos)^2
dtheta=sqrt(tmpn/tmpd)

resol=fltarr(381,2)

for i=0,380 do begin
    en=i+1
    resol(i,0)=en          ;the energy array
    wavel=sqrt(81.805/en)
    stheta=wavel/(2.*dmono)
    if (abs(stheta) lt 1) then begin
       cot_theta=sqrt(1-stheta*stheta)/stheta
    endif else  begin
       cot_theta=0.0
    endelse
    de=2.*en*cot_theta*dtheta
    resol(i,1)=sqrt(de*de+1.2*1.2)
endfor

resfn=resol*8.06544      ; in wavenumbers! (E,Res)
; END RESOLUTION
;_______________________________________________________________________________


spectrum=*self.eigenyptr*scalefactor
eigeny=*self.eigenyptr
freq=*self.freqptr

;need to get the discrete resolutions at each freq
res=fltarr(n_elements(freq))
  for i=0,n_elements(freq)-1 do begin
    P = WHERE((resfn[*,0]) GE freq[i])
    if freq[i] le 2400 then begin
       res[i] =resfn[P[0],1]
    endif else begin
       res[i] =resfn[380,1]
    endelse
  endfor

n1 = n_elements(freqspec)
n2 = n_elements(eigeny)
un1 = 1+bytarr(n1)
un2 = 1+bytarr(n2)
bgarr=freqspec
bgarr=freqspec*B1+B0
sumf = total((un1#eigeny)*exp(-(freqspec#un1-un1#Freq)^2/(un1#res)^2),2)
fundamental = scalefactor*sumf*exp(-(dwf^2)*freqspec)+bgarr
;self.fundamental=fundamental



if (plotover eq 1) then begin
    ; calculate all possible combinations as being
    ; the  sum of frequency and mulitply intensities.
    CombNum = (n_elements(freq)^2 + n_elements(freq))/2
    CombFreq=fltarr(CombNum)
    CombInty=fltarr(CombNum)
       i = 0
        For j = 0,n_elements(freq)-1 do begin
            CombFreq(i) = Freq(j)*2
            CombInty(i) = eigeny(j)^2
            i = i + 1
            For k = (j + 1),n_elements(freq)-1 do begin
                CombFreq(i) = Freq(j) + Freq(k)
                CombInty(i) = eigeny(j) * eigeny(k)
                i = i + 1
            endfor
        endfor


    ;need to get discrete resolutions at each freq
    rescomb=fltarr(n_elements(combfreq))
    for i=0,n_elements(combfreq)-1 do begin
       P = WHERE((resfn[*,0]) GE combfreq[i])
       if combfreq[i] le 2400 then begin
         rescomb[i] =resfn[P[0],1]
       endif else begin
         rescomb[i] =resfn[380,1]
       endelse
    endfor

   n1 = n_elements(freqspec)
   n2 = n_elements(CombInty)

   overtone1=fltarr(n1)
   sum=fltarr(n1)
   un1 = 1+bytarr(n1)
   un2 = 1+bytarr(n2)
   sum = total((un1#CombInty)*exp(-(freqspec#un1-un1#combFreq)^2/(un1#rescomb)^2),2)
   overtone1 = oversf*sum*exp(-(dwf^2)*freqspec)
endif


if plotover then begin
   x=freqspec
   y=overtone1
   data=[[freqspec],[fundamental],[freqspec],[overtone1]]
endif else begin
   data=[[freqspec],[fundamental]]
endelse

void = check_math()
;
;Restore original math error behavior.
;
!except = orig_except

return,data
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function g98log::eigenvec
;from the current mass, vector displacement and eigenvalue
;need to get rid of mass-weighted coordinates and
;calculate the real neutron intensity
mass=(*self.atomnewptr).mass
ixs=(*self.atomnewptr).ixsection
fullvector=*self.vectorptr
eigvec=fltarr(self.dof,self.numberatoms) ; ARRAY SIZE FROM READ
evecs=fltarr(self.dof,self.numberatoms,3)
totty=fltarr(self.dof)
i=0
for i=0,self.numberatoms-1 do begin
   for j=0,self.dof-1 do begin
      eigvec(j,i)=sqrt(fullvector(j,i,0)^2+fullvector(j,i,1)^2+fullvector(j,i,2)^2)
      evecs(j,i,0)=fullvector(j,i,0)*mass(i)
      evecs(j,i,1)=fullvector(j,i,1)*mass(i)
      evecs(j,i,2)=fullvector(j,i,2)*mass(i)
    endfor
endfor

; Renormalize those eigenvectors:
for k=0,self.dof-1 do begin
    tot=0
    esum=0.0
    for j=0,self.numberatoms-1 do begin
       for i=0,2 do begin
         esum=esum+evecs(k,j,i)^2
       endfor

    endfor
    esum=sqrt(esum)
    for j=0,self.numberatoms-1 do begin
       for i=0,2 do begin
         evecs(k,j,i)=evecs(k,j,i)/esum
       endfor
       eigvec(k,j)=sqrt(evecs(k,j,0)^2+evecs(k,j,1)^2+evecs(k,j,2)^2)
       fin=(eigvec(k,j)^2)*ixs(j)*exp(-1.0/mass(j))/mass(j)
       tot=tot+fin
    endfor
    totty[k]=tot
endfor

*self.eigenyptr=totty
*self.massevecsptr=evecs
;XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
;*self.spectrumptr=totty*self.scalefactor
;XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function g98log::loaddata,filename = filename,group_leader = group_leader
if n_elements(group_leader) eq 0 then group_leader = 0L

;catch,errno
;if (errno ne 0) then begin
;  catch,/cancel
;  void = dialog_message(/error,["Problem with your file?",!err_string])
;  return,-1
;endif

if n_elements(filename) eq 0 then begin
  filename = DIALOG_PICKFILE(title = 'Select a G98/G03 output file',$
                             /read,filter = ['*.log','*.out'],$
                             path = self.dataDir, $
                             dialog_parent = group_leader)
endif
if ((file_test(filename) eq 0) or (filename eq '')) then return,-1
self.filename = filename
self.datadir=file_dirname(filename)

   nbig=50000l        ;max size file and text array
      dumstr=''          ;dummy read in variable
      Numlines=0         ;initialize the total number of lines
      startline=0      ;initialize the line we want to read from
      jobline=0          ;this is to extract the G98 job types
      dofline=0          ;this is for the VIB. DEGREES OF FREEDOM
      rottline=0         ;this is for the ROTATIONAL TEMP(K)
; NB 11.6045 K =1meV and  1meV = 8.0655 cm-1, and 1 meV=0.2417 THz
      lowfreqline=0        ;this is for the low frequencies..
      convergeline=0         ;this is for the RMS factors..
      atomline=0         ;this is for the atom counting.
      ORIENTLINE=0       ;Input orientation counting
rotconstants=fltarr(3)   ; hold the rotational temperatures x,y,z in Kelvin
lowfrequencies=fltarr(9)    ; the low frequency vibrations  (floats)
job=''       ; job type array (strings)
dof=0          ; degrees of freedom  (int)
maxforce=fltarr(2)   ; max force obtained, and threshold.
maxforceflag=''        ; did it meet the criteria
rmsforce=fltarr(2)   ; RMS force obtained, and threshold.
rmsforceflag=''        ; did it meet the criteria
maxdisp=fltarr(2)     ; max disp obtained, and threshold.
maxdispflag=''       ; did it meet the criteria
rmsdisp=fltarr(2)     ; RMS disp obtained, and threshold.
rmsdispflag=''       ; did it meet the criteria
numberatoms=0      ; how many atoms in the calculation


      textarr=strarr(nbig)   ;large text array

        openr,lun,filename,/get_lun, error=err
        if (err ne 0) then msg='file isnt there'
        if (err ne 0) then return,msg

        tmp=0
        WHILE (NOT EOF(lun)) DO BEGIN
            readf,lun,dumstr
       textarr(tmp)=dumstr
            if (strcmp(strmid(dumstr,1,2),'# ')) then jobline=tmp
            if (strcmp(strmid(dumstr,1,15),'Deg. of freedom')) then dofline=tmp
            if (strcmp(strmid(dumstr,1,23),'ROTATIONAL TEMPERATURES'))$
            or (strcmp(strmid(dumstr,1,23),'Rotational temperatures'))$
           then rottline=tmp
            if (strcmp(strmid(dumstr,1,15),'Low frequencies')) then lowfreqline=tmp
       if (strcmp(strmid(dumstr,1,13),'Maximum Force')) $
         then convergeline=tmp
       if (strcmp(strmid(dumstr,1,20),'Total atomic charges')) $
          or (strcmp(strmid(dumstr,1,23),'Mulliken atomic charges')) $
         then atomline=tmp
        if (stregex(dumstr,$
                     'and normal coordinates',$
                      /boolean) eq 1) then startline=tmp
         if (stregex(dumstr,$
                     'Input orientation',$
                      /boolean) eq 1) then ORIENTLINE=tmp

       tmp=tmp+1
        ENDWHILE
        numlines=tmp
        free_lun,lun

;extract those things that may prove useful:
if (jobline ne 0) then begin
;    if (n_elements(textarr(jobline)) ge 1) then $
       job=textarr(jobline)
endif

if (dofline ne 0) then begin
    if (n_elements(textarr(dofline)) ge 1) then $
       dofa=(strsplit(textarr(dofline),' ',/extract))
       dof=fix(dofa[3])          ;convert deg. freedom to int.
endif

;dof=9    ; XXXXXXXXXXXXXXXXXXXXXXX


if (dofline ne 0) then begin
    if (n_elements(textarr(rottline)) ge 1) then $
       rota=(strsplit(textarr(rottline),' ',/extract))
    if (n_elements((rota)) ge 4) then $
       rotconstants[0]=rota(3)
    if (n_elements((rota)) ge 5) then $
       rotconstants[1]=rota(4)
    if (n_elements((rota)) ge 6) then $
       rotconstants[2]=rota(5)
endif

if (lowfreqline ne 0) then begin
    if (n_elements(textarr(lowfreqline-1)) ge 1) then $
       lfl=(strsplit(textarr(lowfreqline-1),' ',/extract))
    if (n_elements(lfl) ge 4) then $
       lowfrequencies[0]=lfl(3)
    if (n_elements(lfl) ge 5) then $
       lowfrequencies[1]=lfl(4)
    if (n_elements(lfl) ge 6) then $
       lowfrequencies[2]=lfl(5)
    if (n_elements(lfl) ge 7) then $
       lowfrequencies[3]=lfl(6)
    if (n_elements(lfl) ge 8) then $
       lowfrequencies[4]=lfl(7)
    if (n_elements(lfl) ge 9) then $
       lowfrequencies[5]=lfl(8)

    if (n_elements(textarr(lowfreqline)) ge 1) then $
       lfl=(strsplit(textarr(lowfreqline),' ',/extract))
    if (n_elements(lfl) ge 4) then $
       lowfrequencies[6]=lfl(3)
    if (n_elements(lfl) ge 5) then $
       lowfrequencies[7]=lfl(4)
    if (n_elements(lfl) ge 6) then $
       lowfrequencies[8]=lfl(5)
endif

if (convergeline ne 0) then begin
    maxf=(strsplit(textarr(convergeline),' ',/extract))
    if (n_elements(maxf) ge 2) then $
       maxforce[0]=maxf(2)
    if (n_elements(maxf) ge 3) then $
       maxforce[1]=maxf(3)
    if (n_elements(maxf) ge 4) then $
       maxforceflag=maxf(4)

    rmsf=(strsplit(textarr(convergeline+1),' ',/extract))
    if (n_elements(rmsf) ge 2) then $
       rmsforce[0]=rmsf(2)
    if (n_elements(rmsf) ge 3) then $
       rmsforce[1]=rmsf(3)
    if (n_elements(rmsf) ge 4) then $
       rmsforceflag=rmsf(4)

    maxd=(strsplit(textarr(convergeline+2),' ',/extract))
    if (n_elements(maxd) ge 2) then $
       maxdisp[0]=maxd(2)
    if (n_elements(maxd) ge 3) then $
       maxdisp[1]=maxd(3)
    if (n_elements(maxd) ge 4) then $
       maxdispflag=maxd(4)

    rmsd=(strsplit(textarr(convergeline+3),' ',/extract))
    if (n_elements(rmsd) ge 2) then $
       rmsdisp[0]=rmsd(2)
    if (n_elements(rmsd) ge 3) then $
       rmsdisp[1]=rmsd(3)
    if (n_elements(rmsd) ge 4) then $
       rmsdispflag=rmsd(4)


       if (STRMATCH(rmsdispflag , 'NO' )) then $
          print, 'OOPS!!! best check the RMS DISP. convergence criteria.'
       if (STRMATCH(maxdispflag , 'NO' )) then $
          print, 'OOPS!!! best check the MAX DISP. convergence criteria.'
       if (STRMATCH(rmsforceflag , 'NO' )) then $
          print, 'OOPS!!! best check the RMS FORCE convergence criteria.'
       if (STRMATCH(maxforceflag , 'NO' )) then $
          print, 'OOPS!!! best check the MAX FORCE convergence criteria.'

endif

if (n_elements(textarr(atomline)) ge 1) AND (atomline ne 0) then begin
       i=0
       while $
       ((strcmp(strmid(textarr(atomline+i+2),1,23),'Sum of Mulliken charges'))$
          ne 1) do begin
       tmp=(strsplit(textarr(atomline+i+2),' ',/extract))
       i=i+1
       endwhile
       numberatoms=i
       if (numberatoms eq 0) then print,'You have no atoms? Really?'
endif else begin
       msg='This program seems to be confused. Something in the header lines'
       return,msg
endelse


atomlist=strarr(numberatoms)
newatomlist=strarr(numberatoms)
atomicnumber=intarr(numberatoms)
newatomicnumber=intarr(numberatoms)
ixsection=fltarr(numberatoms)
newixsection=fltarr(numberatoms)
ixsection=fltarr(numberatoms)
mass=fltarr(numberatoms)
covalent=fltarr(numberatoms)
vdw=fltarr(numberatoms)
initialorientation=fltarr(numberatoms,3)

if (orientline ne 0) then begin
         i=0
         while $
         ((strcmp(strmid(textarr(orientline+i+5),1,23),'-----------------------'))$
          ne 1) do begin
         oline=(strsplit(textarr(orientline+i+5),' ',/extract))
         initialorientation[i,0]=oline[3]
         initialorientation[i,1]=oline[4]
         initialorientation[i,2]=oline[5]
         i=i+1
         endwhile

endif

if (n_elements(textarr(atomline)) ge 1) then begin
       i=0
       while $
       ((strcmp(strmid(textarr(atomline+i+2),1,23),'Sum of Mulliken charges'))$
          ne 1) do begin
       tmp=(strsplit(textarr(atomline+i+2),' ',/extract))
       atomlist[i]=tmp[1]
       i=i+1
       endwhile
       if (n_elements(atomlist) eq 0) then print,'You have no atoms? Really?'
endif


startline=startline+2
i=0
dof=1000
;make an arb number of jumps....

symmode=strarr(dof)
freq=fltarr(dof)
redmass=fltarr(dof)
frcc=fltarr(dof)
iri=fltarr(dof)
raman=fltarr(dof)
depolar=fltarr(dof)
vector=fltarr(dof,numberatoms,3)       ;put the displacements in x,y,z form

while (stregex(textarr(startline),'------',/boolean) eq 0) do begin

    symline=(strsplit(textarr(startline),' ',/extract))
    if (n_elements(symline) eq 3) then begin
       symmode[i]=symline[0]
           symmode[i+1]=symline[1]
       symmode[i+2]=symline[2]
       startline=startline+1
      if (strcmp(strmid(textarr(startline),1,14),'Frequencies --')) then begin
          freqline=(strsplit(textarr(startline),' ',/extract))
          freq[i]=freqline[2]
          freq[i+1]=freqline[3]
          freq[i+2]=freqline[4]
          startline=startline+1
       endif
       if (strcmp(strmid(textarr(startline),1,14),'Red. masses --')) then begin
          redline=(strsplit(textarr(startline),' ',/extract))
          redmass[i]=redline[3]
          redmass[i+1]=redline[4]
          redmass[i+2]=redline[5]
          startline=startline+1
       endif
       if (strcmp(strmid(textarr(startline),1,14),'Frc consts  --')) then begin
          frcline=(strsplit(textarr(startline),' ',/extract))
          frcc[i]=frcline[3]
          frcc[i+1]=frcline[4]
          frcc[i+2]=frcline[5]
          startline=startline+1
       endif
       if (strcmp(strmid(textarr(startline),1,14),'IR Inten    --')) then begin
          irline=(strsplit(textarr(startline),' ',/extract))
          irI[i]=irline[3]
          iri[i+1]=irline[4]
           iri[i+2]=irline[5]
          startline=startline+1
       endif
       if (strcmp(strmid(textarr(startline),1,14),'Raman Activ --')) then begin
          ramanline=(strsplit(textarr(startline),' ',/extract))
          raman[i]=ramanline[3]
          raman[i+1]=ramanline[4]
          raman[i+2]=ramanline[5]
          startline=startline+1
       endif
       if (strcmp(strmid(textarr(startline),1,14),'Depolar     --')) then begin
          depolarline=(strsplit(textarr(startline),' ',/extract))
          depolar[i]=depolarline[2]
          depolar[i+1]=depolarline[3]
          depolar[i+2]=depolarline[4]
          startline=startline+1
        endif
          startline=startline+1

       for j=0,numberatoms-1 do begin
         atdisp=(strsplit(textarr(startline),' ',/extract))
         atomicnumber[j]=atdisp[1]
         vector[i,j,0]=atdisp[2]
         vector[i,j,1]=atdisp[3]
         vector[i,j,2]=atdisp[4]
         vector[i+1,j,0]=atdisp[5]
         vector[i+1,j,1]=atdisp[6]
         vector[i+1,j,2]=atdisp[7]
         vector[i+2,j,0]=atdisp[8]
         vector[i+2,j,1]=atdisp[9]
         vector[i+2,j,2]=atdisp[10]
         startline=startline+1
       endfor

       startline=startline+1
       i=i+3
    endif else begin
        symmode[i]=symline[0]
           symmode[i+1]=symline[1]
       startline=startline+1
                freqline=(strsplit(textarr(startline),' ',/extract))
       freq[i]=freqline[2]
       freq[i+1]=freqline[3]
       startline=startline+1
                redline=(strsplit(textarr(startline),' ',/extract))
       redmass[i]=redline[3]
       redmass[i+1]=redline[4]
       startline=startline+1
                frcline=(strsplit(textarr(startline),' ',/extract))
       frcc[i]=frcline[3]
       frcc[i+1]=frcline[4]
       startline=startline+1
                irline=(strsplit(textarr(startline),' ',/extract))
       irI[i]=irline[3]
       iri[i+1]=irline[4]
       startline=startline+1
                ramanline=(strsplit(textarr(startline),' ',/extract))
       raman[i]=ramanline[3]
       raman[i+1]=ramanline[4]
       startline=startline+1
                depolarline=(strsplit(textarr(startline),' ',/extract))
       depolar[i]=depolarline[2]
       depolar[i+1]=depolarline[3]
       depolar[i+2]=depolarline[4]
       startline=startline+2
       for j=0,numberatoms-1 do begin
         atdisp=(strsplit(textarr(startline),' ',/extract))
         atomicnumber[j]=atdisp[1]
         vector[i,j,0]=atdisp[2]
         vector[i,j,1]=atdisp[3]
         vector[i+1,j,0]=atdisp[5]
         vector[i+1,j,1]=atdisp[6]
         vector[i+2,j,0]=atdisp[8]
         vector[i+2,j,1]=atdisp[9]
         startline=startline+1
       endfor
    ENDELSE

if (n_elements(freq) eq 0) then begin
    msg='Either this program is having problems'
    msg=msg+'.....or you dont have any freq info. Try again'
    return,msg
endif
endwhile
;need to populate the parameters

isotopeinfo=(*self.isotopeinfoptr)
for i=0,numberatoms-1 do begin
    j=0
    repeat begin
       test=strsplit((isotopeinfo[j]),',',/extract)
       j=j+1
    endrep until atomlist[i] eq test[0]
    mass[i]=test[1]
    ixsection[i]=test[4]
    covalent[i]=test[6]/1000.  ; in Angstroms
    vdw[i]=test[5]/1000.       ; in Angstroms
endfor

;how many elements?


max=n_elements(where(freq))-1
*self.initialorientationptr=initialorientation
self.rotconstants=rotconstants
self.lowfrequencies=lowfrequencies
self.job=job
self.dof=max+1
self.maxforce=maxforce
self.maxforceflag=maxforceflag
self.rmsforce=rmsforce
self.rmsforceflag=rmsforceflag
self.maxdisp=maxdisp
self.maxdispflag=maxdispflag
self.rmsdisp=rmsdisp
self.rmsdispflag=rmsdispflag
self.numberatoms=numberatoms

atomorig={atomlist:atomlist,atomicnumber:atomicnumber,$
   ixsection:ixsection,mass:mass,covalent:covalent,vdw:vdw}
atomnew={atomlist:atomlist,atomicnumber:atomicnumber,$
   ixsection:ixsection,mass:mass,covalent:covalent,vdw:vdw}
*self.atomorigptr=atomorig
*self.atomnewptr=atomnew

(*self.symmodeptr)=symmode[0:max]
(*self.origfreqptr)=freq[0:max]
(*self.freqptr)=freq[0:max]
(*self.redmassptr)=redmass[0:max]
(*self.frccptr)=frcc[0:max]
(*self.iriptr)=iri[0:max]
(*self.ramanptr)=raman[0:max]
(*self.depolarptr)=depolar[0:max]
(*self.vectorptr)=vector[0:max,*,*]



self.dataLoaded = 1
ret=self->eigenvec()
return,1
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98log::setproperty,temp=temp
end

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98log::getproperty,numberatoms=numberatoms,$
            initcoords=initcoords,              $
            atomlist=atomlist,                  $
            covalent=covalent,                  $
            vdw=vdw,                            $
            eigen=eigen,                        $
            evector=evector,                      $
            massvec=massvec,                    $
            sym=sym,                            $
            freq=freq


if arg_present(numberatoms) then numberatoms=self.numberatoms
if arg_present(initcoords) then initcoords=*self.initialorientationptr
if arg_present(atomlist) then atomlist=(*self.atomnewptr).atomlist
if arg_present(covalent) then covalent=(*self.atomnewptr).covalent
if arg_present(vdw) then vdw=(*self.atomnewptr).vdw
if arg_present(eigen) then eigen =*self.eigenyptr
if arg_present(massvec) then massvec =*self.massevecsptr
if arg_present(evector) then evector =*self.vectorptr
if arg_present(sym) then sym =*self.symmodeptr
if arg_present(freq) then freq =*self.origfreqptr

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function g98log::dataLoaded
return,self.dataLoaded
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
function g98log::init,filename = filename,dataDirectory = dataDirectory
if n_elements(filename) eq 0 then filename = " "
if n_elements(dataDirectory) eq 0 then dataDirectory = $
    'c:\idlstuff\objects'
self.dataDir = dataDirectory
self.filename = filename
self.dataLoaded = 0

;                                                   Van der Waals, Covalent
;name,weigh,COH Scat, Inc Scatt Cross,Total SCatt Cross,milliA,milliA,Abs Cros,linecontinue
isotopeinfo= [$
'H,1.00797,1.7568,80.26,82.02,1200,230,0.3326',$
'D,2.014102,5.592,2.05,7.64,1200,230,5.19E-4',$
'He,4.0026,1.34,0,1.34,1400,930,0.00747',$
'Li,6.942,0.454,0.92,1.37,1820,680,70.5',$
'Be,9.0122,7.63,0.0018,7.63,1700,350,0.0076',$
'B,10.811,3.54,1.7,5.24,2080,830,767.0',$
'C,12.01115,5.551,0.001,5.551,1950,680,0.0035',$
'N,14.0067,11.01,0.5,11.51,1850,680,1.9',$
'O,15.9994,4.232,8.00E-04,4.232,1700,680,1.9E-4',$
'F,18.9984,4.017,8.00E-04,4.018,1730,640,0.0096',$
'Ne,20.183,2.62,0.008,2.628,1540,1120,0.039',$
'Na,22.9898,1.66,1.62,3.28,2270,970,0.53',$
'Mg,24.312,3.631,0.08,3.71,1730,1100,0.063',$
'Al,26.9815,1.495,0.0082,1.503,2050,1350,0.231',$
'Si,28.086,2.163,0.004,2.167,2100,1200,0.171',$
'P,30.9738,3.307,0.005,3.312,2080,750,0.172',$
'S,32.064,1.0186,0.007,1.026,2000,1020,0.53',$
'Cl,35.453,11.5257,5.3,16.8,1970,990,33.5',$
'Ar,39.948,0.458,0.225,0.683,1880,1570,0.675',$
'K,39.102,1.69,0.27,1.96,2750,1330,2.1',$
'Ca,40.08,2.78,0.05,2.83,1973,990,0.43',$
'Sc,44.956,19,4.5,23.5,1700,1440,27.5',$
'Ti,47.9,1.485,2.87,4.35,1700,1470,6.09',$
'V,50.942,0.0184,5.08,5.1,1700,1330,5.08',$
'Cr,51.996,1.66,1.83,3.49,1700,1350,3.05',$
'Mn,54.938,1.75,0.4,2.15,1700,1350,13.3',$
'Fe,55.847,11.22,0.4,11.62,1700,1340,2.56',$
'Co,58.9332,0.779,4.8,5.6,1700,1330,37.18',$
'Ni,58.71,13.3,5.2,18.5,1630,1500,4.49',$
'Cu,63.54,7.485,0.55,8.03,1400,1520,3.78',$
'Zn,65.37,4.054,0.077,4.131,1390,1450,1.11',$
'Ga,69.72,6.675,0.16,6.83,1870,1220,2.75',$
'Ge,72.59,8.42,0.18,8.6,1700,1170,2.2',$
'As,74.9216,5.44,0.06,5.5,1850,1210,4.5',$
'Se,78.96,7.98,0.32,8.3,1900,1220,11.7',$
'Br,79.909,5.8,0.1,5.9,2100,1210,6.9',$
'Kr,83.8,7.67,0.01,7.68,2020,1910,25.0',$
'Rb,85.47,6.32,0.5,6.8,1700,1470,0.38',$
'Sr,87.62,6.19,0.06,6.25,1700,1120,1.28',$
'Y,88.906,7.55,0.15,7.7,1700,1780,1.28',$
'Zr,91.22,6.44,0.02,6.46,1700,1560,0.185',$
'Nb,92.906,6.253,0.0024,6.255,1700,1480,1.15',$
'Mo,95.94,5.67,0.04,5.71,1700,1470,2.48',$
'Tc,99,5.8,0.5,6.3,1700,1350,20.0',$
'Ru,101.07,6.21,0.4,6.6,1700,1400,2.56',$
'Rh,102.905,4.34,0.3,4.6,1700,1450,144.8',$
'Pd,106.4,4.39,0.093,4.48,1630,1500,6.9',$
'Ag,107.87,4.407,0.58,4.99,1720,1590,63.3',$
'Cd,112.4,3.04,3.46,6.5,1580,1690,2520.0',$
'In,114.82,2.08,0.54,2.62,1930,1630,193.8',$
'Sn,118.69,4.871,0.022,4.892,2170,1460,0.626',$
'Sb,121.75,3.9,0.007,3.9,2200,1460,4.91',$
'Te,127.6,4.23,0.09,4.32,2060,1470,4.7',$
'I,126.9044,3.5,0.31,3.81,2150,1400,6.15',$
'Xe,131.3,2.96,0,0,2160,1980,23.9',$
'Cs,132.905,3.69,0.21,3.9,1700,1670,29.0',$
'Ba,137.34,3.23,0.15,3.38,1700,1340,1.1',$
'La,138.91,8.53,1.13,9.66,1700,1870,8.97',$
'Ce,140.12,2.94,0.001,2.94,1700,1830,0.63',$
'Pr,140.907,2.64,0.015,2.66,1700,1820,11.5',$
'Nd,144.24,7.43,9.2,16.6,1700,1810,50.5',$
'Pm,145,20,1.3,21.3,1700,1800,168.4',$
'Sm,150.35,0.422,39,39,1700,1800,5922.0',$
'Eu,151.96,6.57,2.5,9.2,1700,1990,4530.0',$
'Gd,157.25,29.3,151,180,1700,1790,49700.0',$
'Tb,158.925,6.84,0.004,6.84,1700,1760,23.4',$
'Dy,162.5,35.9,54.4,90.3,1700,1750,994.0',$
'Ho,164.93,8.06,0.36,8.42,1700,1740,64.7',$
'Er,167.26,7.63,1.1,8.7,1700,1730,159.0',$
'Tm,168.934,6.28,0.1,6.38,1700,1720,100.0',$
'Yb,173.04,19.42,4,23.4,1700,1940,34.8',$
'Lu,174.97,6.53,0.7,7.2,1700,1720,74.0',$
'Hf,178.49,7.6,2.6,10.2,1700,1570,104.1',$
'Ta,180.948,6,0.01,6.01,1700,1430,20.6',$
'W,183.85,2.97,1.63,4.6,1700,1370,18.3',$
'Re,186.2,10.6,0.9,11.5,1700,1350,89.7',$
'Os,190.2,14.4,0.3,14.7,1700,1370,16.0',$
'Ir,192.2,14.1,0,14,1700,1320,425.0',$
'Pt,195.09,11.58,0.13,11.71,1720,1500,10.3',$
'Au,196.967,7.32,0.43,7.75,1660,1500,98.65',$
'Hg,200.59,20.24,6.6,26.8,1550,1700,372.3',$
'Tl,204.37,9.678,0.21,9.89,1960,1550,3.43',$
'Pb,207.19,11.115,0.003,11.118,2020,1540,0.171',$
'Bi,208.98,9.148,0.0084,9.156,1700,1540,0.0338',$
'Po,210,0,0,00,1700,1680',$                   ,
'At,214.9986412,0,0,0,1700,1700,0',$
'Rn,222.0175705,0,0,12.6,1700,2400,0',$
'Fr,224,0,0,0,1700,2000,0',$
'Ra,226.0254,13,0,13,1700,1900,12.8',$
'Ac,227,0,0,0,1700,1880,0',$
'Th,2332.038,13.36,0,13.36,1700,1790,7.37',$
'Pa,232.5,10.4,0.1,10.5,1700,1610,200.6',$
'U,238.03,8.903,0.005,8.908,1860,1580,7.57',$
'Np,237,14,0.5,14.5,1700,1550,175.9',$
'Pu,239.0522,0,0,0,1700,1530,0',$
'Am,242,8.7,0.3,9,1700,1510,75.3',$
'Cm,247.5,0,0,0,1700,1500,0']
self.isotopeinfoptr=ptr_new(isotopeinfo)
self.rotconstants=fltarr(3)   ; hold the rotational temperatures x,y,z in Kelvin
self.lowfrequencies=fltarr(9)    ; the low frequency vibrations  (floats)
self.job=''       ; job type array (strings)
self.dof=0          ; degrees of freedom  (int)
self.maxforce=fltarr(2)   ; max force obtained, and threshold.
self.maxforceflag=''        ; did it meet the criteria
self.rmsforce=fltarr(2)   ; RMS force obtained, and threshold.
self.rmsforceflag=''        ; did it meet the criteria
self.maxdisp=fltarr(2)     ; max disp obtained, and threshold.
self.maxdispflag=''       ; did it meet the criteria
self.rmsdisp=fltarr(2)     ; RMS disp obtained, and threshold.
self.rmsdispflag=''       ; did it meet the criteria
self.numberatoms=0      ; how many atoms in the calculation
self.initialorientationptr=ptr_new(/allocate_heap)
self.atomorigptr=ptr_new(/allocate_heap)
self.atomnewptr=ptr_new(/allocate_heap)
self.symmodeptr=ptr_new(/allocate_heap)
self.origfreqptr=ptr_new(/allocate_heap)
self.freqptr=ptr_new(/allocate_heap)
self.redmassptr=ptr_new(/allocate_heap)
self.frccptr=ptr_new(/allocate_heap)
self.iriptr=ptr_new(/allocate_heap)
self.ramanptr=ptr_new(/allocate_heap)
self.depolarptr=ptr_new(/allocate_heap)
self.vectorptr=ptr_new(/allocate_heap)
self.eigenyptr=ptr_new(/allocate_heap)   ; the eigenvector for each mode
self.massevecsptr=ptr_new(/allocate_heap)   ; the mass weighted normalized eigenvector
if self.filename ne " " then begin
 ret =  self->loadData(filename = self.filename)
 return,1
endif

end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98log__define
; Definition module for g98log class
define = {g98log,                   $
          dataDir:"",               $
          dataLoaded:0,             $
          filename:"",              $
          isotopeinfoptr:ptr_new() ,$    ;  CONTAINS ALL ATOMINFO
          rotconstants:fltarr(3),   $    ; hold the rotational temperatures x,y,z in Kelvin
          lowfrequencies:fltarr(9), $    ; the low frequency vibrations  (floats)
          job:'',                   $    ; job type array (strings)
          dof:0,                    $    ; degrees of freedom  (int)
          maxforce:fltarr(2),       $    ; max force obtained, and threshold.
          maxforceflag:'',          $    ; did it meet the criteria
          rmsforce:fltarr(2),       $    ; RMS force obtained, and threshold.
          rmsforceflag:'',          $    ; did it meet the criteria
          maxdisp:fltarr(2),        $    ; max disp obtained, and threshold.
          maxdispflag:'',           $    ; did it meet the criteria
          rmsdisp:fltarr(2),        $    ; RMS disp obtained, and threshold.
          rmsdispflag:'',           $    ; did it meet the criteria
          numberatoms:0 ,           $    ; how many atoms in the calculation
          initialorientationptr:ptr_new(),$
          atomorigptr:ptr_new(),   $
          atomnewptr:ptr_new(),    $
          symmodeptr:ptr_new() ,$
          origfreqptr:ptr_new(),$
          freqptr:ptr_new(),    $
          redmassptr:ptr_new(), $
          frccptr:ptr_new(),    $
          iriptr:ptr_new(),     $
          ramanptr:ptr_new(),   $
          depolarptr:ptr_new(), $
          vectorptr:ptr_new()  ,$
          eigenyptr:ptr_new(),   $  ; the eigenvector for each mode
          spectrumptr:ptr_new(),   $  ;
          massevecsptr:ptr_new()   $  ;mass weighted normalized vector
}

return
end
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
pro g98_ex,g03=g03,g98=g98

if g03 then filename = "c:\idlstuff\objects\inma.out"
if g98 then filename = "c:\idlstuff\objects\pyz.log"

o = obj_new('g98log',filename = filename)
print,o->dataloaded()
help,o
obj_destroy,o
end
