PRO dpoint, pause = pause ;, tof = toff

	; Purturbs the zenith angle of a synthetic spectrum
	; & retrieves the difference in the pbp
	; assumes we're in a working rdrv retrieval directory
		
	FORWARD_FUNCTION readstat, readctl, setfilenames, readspec, mksyn, readt15

	; Set up file names from LAPACK version sfit2 run 
	rc = setfilenames( fn )
		
	; Use default ILS for all calcs
	cmd = 'cp ' + fn.ilsfit + ' ' + fn.ilsnam
	SPAWN, cmd	

	; Use default cinput for all calcs
	cmd = 'cp ' + fn.cindef + ' ' + fn.cinnam
	SPAWN, cmd
	
	; Use default fasc.inp for all calcs
	cmd = 'cp ' + fn.fasdef + ' ' + fn.fasnam
	SPAWN, cmd
	
	; Use default fasc.prf for all calcs
	cmd = 'cp ' + fn.prfdef + ' ' + fn.prfnam
	SPAWN, cmd
	
	; Set default rdrv.ctl at start
	cmd = 'cp ' + fn.rdfdef + ' ' + fn.rdfnam
	SPAWN, cmd

	fpert = 0.2 				; set in "fasc.inp.per"
	junk=''
	
	PRINT, ' The zenith angle perturbation value [deg] : ', fpert
	PRINT, ' Default fasc.inp   : ', fn.fasdef
	PRINT, ' Perturbed fasc.inp : ', fn.fasper
	READ, junk, PROMPT=' Please confirm that fasc.inp.per exists with the perturbation to the sza = sza+0.2. Hit enter to continue or q(uit): '
    IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
  PRINT,' Continuing..'
		

	; Define codes to do the work
  fast = '/project/ya4/bin/fastc.k fasc'
  sfit = '/project/ya4/bin/sfit2.k -i -v 394xo' 
  rdrv = '/project/ya4/sfitsrc/394lp/rdrv39'

	; STEP 1
	
	; Make the default spectrum & do a retrieval
	; Make the dummy spectrum
	mksyn, file = fn.t15one
	
	; Create a forward model rdrv.ctl file from the existing one
	rc = 0
	rc = readctl( ctl, fn.rdfnam )
	
	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15one
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, '.FALSE.'		
	FREE_LUN, lun
	
	; Compute the default vmr, column
	SPAWN, fast	 
	SPAWN, rdrv

	junk=''
	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT=' Finished default spectrum - Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT,' Finished default spectrum..'


	; Compute a retrieval from the default spectrum
	; Make a retrieval rdrv.ctl
	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, 'K.out'		
		PRINTF, lun, '.TRUE.'
		PRINTF, lun, 'Sa.complete'		
	FREE_LUN, lun

	; Compute the default retrieval, vmr, column ...
	SPAWN, rdrv

	;; Save the necessary files
	;cmd = 'cp PRFS.out PRFS.out.def'
	;SPAWN, cmd
	;cmd = 'cp statevec statevec.def'
	;SPAWN, cmd
	;cmd = 'cp cdetail cdetail.def'
	;SPAWN, cmd
  ;cmd = 'cp pbpfile pbpfile.def'
  ;SPAWN, cmd

;read in the unperturbed pbp file
  
  unpertpbp=readspec(fn.pbfile)

	; Read in Statevector file & get neutral vmr
	rc = 0
	rc = readstat( stdef, fn.stfile )
	
	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT=' Finished retrieval of default - Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT, ' Finished retrieval of default...'


	; STEP 2
	
	; Perturb the pointing by re-calculating the airmass with a new viewing angle
	;  and make a synthetic spectrum

	; Make a forward model rdrv.ctl
	PRINT, 'Overwriting : ', fn.rdfnam
	OPENW, lun, fn.rdfnam, /GET_LUN
		PRINTF, lun, ctl.binf
		PRINTF, lun, ctl.ptfl
		PRINTF, lun, ctl.msfl
		PRINTF, lun, ctl.mxfl
		PRINTF, lun, fn.t15one
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'
		PRINTF, lun, '.FALSE.'		
	FREE_LUN, lun

	; Use perturbed fasc.inp
	cmd = 'cp ' + fn.fasper + ' ' + fn.fasnam
	SPAWN, cmd

	; Compute the atmosphere and retrieve new profile
	SPAWN, fast	 

;	des = ''
;	READ, des, PROMPT='  Did you edit cinput.per with zenith angle code from perturbed fasc.ms? (y/n) : '
;	IF( STRUPCASE(des) NE 'Y' ) THEN STOP
	
	 ; Update the sza with zenith angle from perturbed fasc.ms
  openr,1,'fasc.ms'
  readf,1,junk
  close,1

  tempvalue=strsplit(junk,/extract) 
  newsza=fix(tempvalue(0))
  sedcommand='sed s/NNNNN/'+tempvalue(0)+'/g cinput.def.NNNNN > cinput'
  spawn, sedcommand
  
	
	SPAWN, rdrv

	IF( KEYWORD_SET( pause )) THEN BEGIN
		READ, junk, PROMPT='Finished perturbed spectrum -  Hit enter to continue or q(uit): '
		IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
	ENDIF ELSE PRINT, 'Finished perturbed spectrum...'

;	; Retrieve vmr/column with original pointing forward model from perturbed spectrum
	
;	; Make a retrieval rdrv.ctl	
;	PRINT, 'Overwriting : ', fn.rdfnam
;	OPENW, lun, fn.rdfnam, /GET_LUN
;		PRINTF, lun, ctl.binf
;		PRINTF, lun, ctl.ptfl
;		PRINTF, lun, ctl.msfl
;		PRINTF, lun, ctl.mxfl
;		PRINTF, lun, fn.t15syn
;		PRINTF, lun, '.TRUE.'
;		PRINTF, lun, '.TRUE.'
;		PRINTF, lun, 'K.out'		
;		PRINTF, lun, '.TRUE.'
;		PRINTF, lun, 'Sa.complete'		
;	FREE_LUN, lun
 
	; Use default fasc.inp
;	cmd = 'cp ' + fn.fasdef + ' ' + fn.fasnam
;	SPAWN, cmd

	; Use default cinput
;	cmd = 'cp ' + fn.cindef + ' ' + fn.cinnam
;	SPAWN, cmd

;	SPAWN, fast
;	SPAWN, rdrv

; ; Read in Statevector file 
; rc = 0
; rc = readstat( stat, fn.stfile )

; IF( rc NE 0 ) THEN BEGIN
;   PRINTF, -2,'Could not read st file: ', fn.stfile
;   STOP
; ENDIF

 ;read in pbp
	newpbp=readspec(fn.pbfile)  
	

;calculate Kb : change in spectrum/change in sza (old-new)/perturbation

dif = unpertpbp - newpbp

;set up the Kb matrix
 temp=size(unpertpbp)
 npoints=temp(1)
 Kb=fltarr(npoints)

Kb=dif/fpert

;print out the file. This will be used in errorcalcncar.
openw,1,'Kb_pointing.cov'
for i=0, npoints-1 do printf,1,Kb(i)
close,1



END