;
; computes volumne ratio of Si and C for given porosity Vvac matching 3.6 < [Si]/[C] < 4.
; enter first Vvac = and then later vsi and vc using sugested output values:
;
; run DarVol.pro
;
; Result:
; [Si]/[C] Vsi    Vc     por     Veis
;    3.6   0.59   0.41   0.00    0.
;    3.6   0.57   0.38   0.05    0.
;    3.7   0.53   0.37   0.10    0.
;    3.6   0.48   0.32   0.20    0.
;    3.7   0.41   0.29   0.30    0.
;    3.6   0.36   0.24   0.40    0.
;    3.6   0.30   0.20   0.50    0.
; -------------------------------------------------------------
;

        rhsi   = 3.4
	wmolsi = 134.5
	rhc    = 1.6
	wmolc  = 12.
	
        Vvac = 0.05


	Vsi  = findgen (10)/10. * 0.7 + 0.3
        Vc   = 1. - vvac - Vsi
        Vd   = Vc + VSi + Vvac
        rhsi    = rhsi/(1.+Vvac)
        rhc     = rhc /(1.+Vvac)
;     mass in 1 gramm of dust in aC and Si sub-particles;
        md_ac     = (Vc*rhc)   /(Vc*rhc + Vsi*rhsi)
        md_Si     = (Vsi*rhsi) /(Vc*rhc + Vsi*rhsi)
        md_aCSi   = md_aC / md_Si
        abud_acSi = md_aCSi *wmolsi/wmolc


del = abs(abud_acSi -3.6)
ii  = (where(del  eq min(del))) (0)
print, ' abud_acSi(ii,ii-1)', abud_acSi(ii), abud_acSi(ii-1)
	if  max(vd) ne min(vd)  and max(vd) ne 1 then stop, '** Check : total volumne ne 1', vd
	if  md_ac(ii)+md_si(ii)             ne 1 then stop, '** Check : totalmass not 1g', md_ac(ii)+md_si(ii)

print, format='(a30, f7.1, 3f7.2)', ' *** Volumen: [Si]/[C], Vsi, Vc, por = ', abud_acSi(ii), Vsi(ii), Vc(ii), Vvac
print, format='(a30, f7.1, 3f7.2)', ' *** Volumen: [Si]/[C], Vsi, Vc, por = ', abud_acSi(ii-1), Vsi(ii-1), Vc(ii-1), Vvac
;
; -------------------
;

 print, ' select by hands : Vsi, Vc = ?'

   vsi = 0.55
   vc  = 0.4
 
        Vd   = Vc + VSi + Vvac
	if  max(vd) ne 1. then print, '** Check : total volumne ne 1', vd

           rhsi    = rhsi/(1.+Vvac)
           rhc     = rhc /(1.+Vvac)
	   

c     mass in 1 gramm of dust in aC and Si sub-particles;
        md_ac     = (Vc*rhc)   /(Vc*rhc + Vsi*rhsi)
        md_Si     = (Vsi*rhsi) /(Vc*rhc + Vsi*rhsi)
        md_aCSi   = md_aC / md_Si
        abud_acSi = md_aCSi *wmolsi/wmolc
	print, format='(a30, f7.1, 3f7.2)', ' *** Final volumen: [Si]/[C], Vsi, Vc, por = ', abud_acSi, Vsi, Vc, Vvac

end