;This Script will open up a 'test.qdp' file output from Andy Rasmussen's 
;vis_bi_ccd.c code and pull out the following quantities
;
;Z  =  lab g1 CCD z (from frontside) [\gmm]
;D  =  lab g2 doping density [cm\u-3\d]
;E  =  lab g3 field strength [kV cm\u-1\d]
;V  =  lab g4 field potential [V]
;R  =  lab g5 diffusion scale [\gmm]

pro vis_ccd_profile,filename=filename

If (not Keyword_Set(filename)) then $
Filename = '/nfs/slac/g/ki/ki08/lsst/lsstgeant4/050117_lsst_vis_bi_ccd/test.qdp'
If (not Keyword_Set(Dd)) then Dd = 200.e-6
If (not Keyword_Set(TvFlag)) then TvFlag=0
PSFileName='/nfs/slac/g/ki/ki08/lsst/lsstgeant4/plots/DiffusionRadius.ps'
tvflag=1
;Constants:  
k=double(1.3806e-23) ;cm^2 kg / s^2 K
q=double(1.609e-19)  ;Coulombs 
T=40.	 	     ;Kelvin
Vb=30.		     ;Volts
PixelSize=10.	     ;Microns
;readcol,FileName,Z,D,E,V,R,format='f,f,f,f,f'
;Fitting Parameter

s=5.1
;Get Diffusion Radius from G.G Pavlov and J.A. Nouse Nucl. Instr. and Meth. 
;A 428 (1999) 348
Zd=Dd*(findgen(10000)/10000)
Rd=(Zd*Dd*2*k*T/(q*Vb))^.5*(s+(Dd/Zd)*alog(Dd/(Dd-Zd)))^.5
RdMicrons=Rd*10.^6
RdMicronsMax=15
RdMicronsMin=0

ZdivDd=Zd/Max(Zd)

if tvflag eq 1 then begin
  set_plot,'x'
endif else begin
   Set_plot,'ps'
   device,filename=PSFileName
   loadct,39                       ;load color table 39
   device,/color                   ;allow color on the postscript
   device,ysize=7.0,/inches        ;Height of plot in y
   device,xsize=8.0,/inches        ;Width of plot in x
   device,yoffset=0.0,/inches      ;Y position of lower left corner
   device,xoffset=0.0,/inches

   !p.charsize=1.5
   !p.thick=4.5
   !p.charthick=4
   !x.thick=4
   !y.thick=4
   white='FFFFFF'x
   black='000000'x
   red='FF0000'x
   green='00FF00'x
   blue='0000FF'x
endelse

plot,ZdivDd,RdMicrons,xtitle='z/d!ld!N',ytitle='Microns',position=[0.1,0.1,0.9,0.9],/norm,$
     yrange=[RdMicronsMin,RdMicronsMax]
oplot,ZdivDd,Reverse(R),linestyle=3
axis,yaxis=1,ytitle='Pixels',$
     yrange=[RdMicronsMin/PixelSize,RdMicronsMax/PixelSize]
if tvflag ne 1 then begin
  device,/close
  set_plot,'x'
endif
stop 

end