How can I deconvolve a spectrum to get the information about each of the deconvolved spectra?

I have a spectrum with two peaks. I want to deconvolve it by gaussian so that I get the parametres (e.g. peak position, area, FWHM, Height) of each deconvolved spectrum.
I have tried the following:

data={{400, 0.04}, {401, 0.039}, {402, 0.04}, {403, 0.039}, {404, 
  0.039}, {405, 0.039}, {406, 0.04}, {407, 0.04}, {408, 0.04}, {409, 
  0.04}, {410, 0.04}, {411, 0.04}, {412, 0.04}, {413, 0.041}, {414, 
  0.041}, {415, 0.042}, {416, 0.041}, {417, 0.041}, {418, 
  0.041}, {419, 0.042}, {420, 0.042}, {421, 0.042}, {422, 
  0.042}, {423, 0.042}, {424, 0.043}, {425, 0.043}, {426, 
  0.043}, {427, 0.044}, {428, 0.045}, {429, 0.046}, {430, 
  0.047}, {431, 0.048}, {432, 0.049}, {433, 0.049}, {434, 
  0.051}, {435, 0.052}, {436, 0.053}, {437, 0.055}, {438, 
  0.057}, {439, 0.059}, {440, 0.061}, {441, 0.063}, {442, 
  0.066}, {443, 0.068}, {444, 0.071}, {445, 0.074}, {446, 
  0.076}, {447, 0.08}, {448, 0.083}, {449, 0.088}, {450, 0.092}, {451,
   0.097}, {452, 0.101}, {453, 0.105}, {454, 0.11}, {455, 
  0.115}, {456, 0.12}, {457, 0.125}, {458, 0.131}, {459, 0.136}, {460,
   0.142}, {461, 0.148}, {462, 0.154}, {463, 0.161}, {464, 
  0.168}, {465, 0.175}, {466, 0.183}, {467, 0.192}, {468, 
  0.201}, {469, 0.211}, {470, 0.221}, {471, 0.233}, {472, 
  0.247}, {473, 0.26}, {474, 0.276}, {475, 0.294}, {476, 0.311}, {477,
   0.329}, {478, 0.349}, {479, 0.369}, {480, 0.391}, {481, 
  0.413}, {482, 0.437}, {483, 0.46}, {484, 0.484}, {485, 0.508}, {486,
   0.531}, {487, 0.554}, {488, 0.578}, {489, 0.6}, {490, 0.622}, {491,
   0.643}, {492, 0.663}, {493, 0.681}, {494, 0.699}, {495, 
  0.715}, {496, 0.731}, {497, 0.745}, {498, 0.759}, {499, 
  0.773}, {500, 0.786}, {501, 0.8}, {502, 0.815}, {503, 0.832}, {504, 
  0.85}, {505, 0.87}, {506, 0.894}, {507, 0.921}, {508, 0.95}, {509, 
  0.982}, {510, 1.019}, {511, 1.059}, {512, 1.102}, {513, 
  1.149}, {514, 1.199}, {515, 1.25}, {516, 1.302}, {517, 1.356}, {518,
   1.409}, {519, 1.461}, {520, 1.511}, {521, 1.557}, {522, 
  1.597}, {523, 1.63}, {524, 1.656}, {525, 1.674}, {526, 1.683}, {527,
   1.683}, {528, 1.672}, {529, 1.651}, {530, 1.62}, {531, 
  1.579}, {532, 1.528}, {533, 1.469}, {534, 1.403}, {535, 1.33}, {536,
   1.25}, {537, 1.165}, {538, 1.08}, {539, 0.996}, {540, 0.913}, {541,
   0.83}, {542, 0.752}, {543, 0.677}, {544, 0.606}, {545, 
  0.541}, {546, 0.481}, {547, 0.426}, {548, 0.375}, {549, 0.33}, {550,
   0.289}, {551, 0.252}, {552, 0.219}, {553, 0.19}, {554, 
  0.165}, {555, 0.143}, {556, 0.123}, {557, 0.107}, {558, 
  0.092}, {559, 0.078}, {560, 0.067}, {561, 0.058}, {562, 0.05}, {563,
   0.042}, {564, 0.036}, {565, 0.031}, {566, 0.026}, {567, 
  0.024}, {568, 0.02}, {569, 0.018}, {570, 0.015}, {571, 0.012}, {572,
   0.011}, {573, 0.01}, {574, 0.008}, {575, 0.007}, {576, 
  0.006}, {577, 0.005}, {578, 0.004}, {579, 0.003}, {580, 
  0.003}, {581, 0.003}, {582, 0.002}, {583, 0.001}, {584, 
  0.001}, {585, 0.001}, {586, 0.001}, {587, 0}, {588, 
  0}, {589, -0.001}, {590, -0.001}, {591, -0.001}, {592, -0.001}, 
{593, -0.001}, {594, -0.001}, {595, -0.002}, {596, -0.002}, {597, 
-0.002}, {598, -0.002}, {599, -0.002}, {600, -0.002}, {601, -0.002}, 
{602, -0.002}, {603, -0.003}, {604, -0.003}, {605, -0.003}, {606, 
-0.003}, {607, -0.003}, {608, -0.003}, {609, -0.003}, {610, -0.003}, 
{611, -0.003}, {612, -0.003}, {613, -0.003}, {614, -0.003}, {615, 
-0.003}, {616, -0.002}, {617, -0.003}, {618, -0.003}, {619, -0.003}, 
{620, -0.003}, {621, -0.003}, {622, -0.003}, {623, -0.003}, {624, 
-0.003}, {625, -0.003}, {626, -0.003}, {627, -0.003}, {628, -0.003}, 
{629, -0.003}, {630, -0.003}, {631, -0.003}, {632, -0.003}, {633, 
-0.003}, {634, -0.003}, {635, -0.003}, {636, -0.003}, {637, -0.003}, 
{638, -0.003}, {639, -0.003}, {640, -0.003}, {641, -0.003}, {642, 
-0.003}, {643, -0.003}, {644, -0.003}, {645, -0.003}, {646, -0.003}, 
{647, -0.003}, {648, -0.003}, {649, -0.003}, {650, -0.003}}

fit = 

    NonlinearModelFit[
      data, 
      y0 + A1/(w1*[Sqrt](2 [Pi]))*Exp[-1/2*((t - xc1)/w1)^2] + 
       A2/(w2*[Sqrt](2 [Pi]))*Exp[-1/2*((t - xc2)/w2)^2], {{y0, 
        0.1}, {w1, 50}, {xc1, 493}, {w2, 20}, {xc2, 526}, {A1, 30}, {A2, 
        48}}, t]

    fitvalue = fit["BestFitParameters"]
    y0 = y0 /. fit["BestFitParameters"][[1]];
    xc1 = xc1 /. fit["BestFitParameters"][[3]];
    w1 = w1 /. fit["BestFitParameters"][[2]];
    xc2 = xc2 /. fit["BestFitParameters"][[5]];
    w2 = w2 /. fit["BestFitParameters"][[4]];
    A1 = A1 /. fit["BestFitParameters"][[6]];
    A2 = A2 /. fit["BestFitParameters"][[7]];
    
    peak1 = N[
       Table[{t, 
         y0 + A1/(w1*[Sqrt](2 [Pi]))*Exp[-1/2*((t - xc1)/w1)^2]}, {t, 
         400,650,1}]];
    peak2 = N[
       Table[{t, 
         y0 + A2/(w2*[Sqrt](2 [Pi]))*Exp[-1/2*((t - xc2)/w2)^2]}, {t, 
         400,650, 1}]];
    sumpeak = 
      N[Table[{t, fit[t]}, {t, 400,650, 1}]];
    Show[ListPlot[data, 
      PlotStyle -> {Blue, Thin}], 
     ListLinePlot[sumpeak, PlotStyle -> {Red, Thick}], 
     ListLinePlot[peak1, PlotStyle -> {Green, Thin}, 
      PlotRange -> {{400, 650}, {0, 2}}], 
     ListLinePlot[peak2, PlotStyle -> {Green, Thin}], 
     PlotRange -> {{400, 650}, {0, 2}}]

(where A1 and A2 are the area of the deconvoluted spectra an w1, w2 are the FWHM and xc1, xc2 are the peak positions of corresponding spectra)
Now, sometime it works, but not always. When the guess value are very close to actual fitted one, then it works, otherwise it does not.
So, how can I modify this or how to solve the problem?

Thank you in advance !