PRO dlshap, pause = pause ;, tof = toff

   ; edit cinput ieap = 1 and iephs=0 call it cinput.ils


	; Substitutes a 20% loss in modulation by replacing a tabular file of
	; mod eff values
	; No phase
	
	; does not fit mod eff
		
	; Runs one version of rdev39 to make a spectrum with a preturbed parameter
	; Then does a retrieval using the default version
	
	; Requires a valid t15 file...
	
	
	FORWARD_FUNCTION readstat, readctl, setfilenames, mksyn, readt15

	; Set up file names from LAPACK version sfit2 run 
	rc = setfilenames( fn )

	uflag = 0
	ugas = 0
	uval = 0.0
	junk=''
	
	READ, junk, PROMPT=' This program needs cinput.ils, where IEAP has been changed to 1 and IEPHS = 0. Please confirm that cinput.ils exists. Hit enter to continue or q(uit): '
    IF( STRUPCASE( junk ) EQ 'Q' )THEN STOP, '   ...Quit by command.'
  PRINT,' Continuing..'
  
		
	; Use 100% ILS for default
	; not linefit
	cmd = 'cp ' + fn.ils100 + ' ' + fn.ilsnam
	SPAWN, cmd	

	; Use default cinput for all calcs
	cmd = 'cp ' + fn.cinils + ' ' + 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

	; 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
	
	; Do a retrieval on a default synthetic spectrum and save the profile & column
	
	; Create a forward model rdrv.ctl file from the existing default one
	rc = 0
	rc = readctl( ctl, fn.rdfnam )

	; Make the forward model rdrv.ctl
	mksyn, file = fn.t15one
	
	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 airmass & spectrum
	SPAWN, fast
	SPAWN, rdrv

	junk=''
	IF( KEYWORD_SET( pause )) THEN BEGIN
		PRINT,''
		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...continue'
	
	; 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
	
	; Read in Statevector file & get neutral vmr
	rc = 0
	rc = readstat( stdef, fn.stfile )
	
	;read in the unperturbed pbp file
  
  unpertpbp=readspec(fn.pbfile)

	
	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 
	; Make a perturbed spectrum & do a retrieval
	
	cmd = 'cp ' + fn.ilsper + ' ' + fn.ilsnam
	SPAWN, cmd

	; Make the forward model rdrv.ctl - overwrite t15syn W/ perturbation parameters	
	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.'		; retrieval flag
		PRINTF, lun, fn.t15syn
		PRINTF, lun, '.FALSE.'		; write k-mat flag
		PRINTF, lun, '.FALSE.'		; write Sa flag
	;	PRINTF, lun, uflag, ugas, uval
	FREE_LUN, lun
		
	; Make perturbed synthetic spectrum	
	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...'
	
	;; 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 ILS	
	;cmd = 'cp ' + fn.ils100 + ' ' + fn.ilsnam
	;SPAWN, cmd

	; Do default retrieval on perturbed spectrum
	;SPAWN, rdrv

	; Read in Statevector file & get purturbed vmr
	;rc = 0
	;rc = readstat( stat, fn.stfile )
	
	;read in pbp
  newpbp=readspec(fn.pbfile)  
  
  ;calculate Kb : change in spectrum/change in Temperature (old-new)/perturbation

dif = unpertpbp - newpbp
fpert=20 ;20% perturbation in the line shape

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

Kb=dif/fpert

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


END