Null(missing) wave during a fitting for user-defined function

I cannot figure out why this function, Do_fit(), does not works. It makes "while executing a wave assignment, .... on a null (missing) wave"

coefwave = coefficient wave
I_exp = experimental data to be fitted
q_x = experimental X wave
I_fit = y wave destination

Function Do_Fit()
    wave coefwave, I_exp, q_x
    Duplicate/O I_exp I_fit
    I_fit = 0
   
    FuncFit/NTHR=0 Fitting_Engine coefwave I_exp /X=q_x /I=1 /D= i_fit

End
Possibly I_fit is not being created because I_exp does not exist.

To test this, add this statement after the wave declarations:
if (!WaveExists(I_exp))
  Abort "ERROR"
endif


I also recommend that you, and anyone else doing Igor programming, turn on "Debug on Error". For details, execute this:
DisplayHelpTopic "Debugging on Error"


In this case, the "NVAR SVAR WAVE Checking" feature will break into the debugger if coefwave, I_exp, or q_x do not exist.

Sometimes when there is an error and Igor breaks into the debugger, the arrow will point to the line AFTER the line that caused the error, so you need to be aware of that possibility.
Thank you, there was no problem for the recall of waves but seems fitting_engine has a problem.

After  return Fi_tot , the fitting process does not works and no error massages appears in the Debug mode.

The governing equation to fit is Fi_tot.

Function Fitting_Engine(pw,xw)
    wave pw
    variable xw


// PERFECT UNIT CELL //
    wave FF_uc, FF_uc_Y, FF_uc_Zr, FF_uc_O
    nvar uc_occ_Y, uc_occ_Zr,  uc_occ_O, ci, u_Y, u_Zr, u_O
    wave ASF_Y, ASF_Zr, ASF_O, uc_position_Y, uc_position_Zr, uc_position_O
   
    FF_uc_Y = atomic_frac_Y * ASF_Y * uc_occ_Zr * exp(ci*xw*(uc_position_Y[0]))*exp(-0.5*(xw*u_Y)^2)
    FF_uc_Zr = atomic_frac_Zr * ASF_Zr * uc_occ_Zr * exp(ci*xw*(uc_position_Zr[0]))*exp(-0.5*(xw*u_Zr)^2)
    FF_uc_O = 2*atomic_frac_O * ASF_O * uc_occ_O * exp(ci*xw*(uc_position_O[0]))*exp(-0.5*(xw*u_O)^2)

        FF_uc = FF_uc_Y + FF_uc_Zr + FF_uc_O   

// CRYSTAL TRUNCATION ROD //
    wave FF_ctr
   
        FF_ctr = 1 / (1-exp(-ci*xw*c))

// ADSORBED ATOMS //
    variable i_duc
    wave duc_position_O
    Duplicate/O duc_position_O temp_duc_position_O
    Duplicate/O duc_position_O temp_duc_position_Y
    Duplicate/O duc_position_O temp_duc_position_Zr
       
    temp_duc_position_O[0] = (1 + pw[4]) * c
    temp_duc_position_Y[0] = (1 + pw[4]) * (c/2)
    temp_duc_position_Zr[0] = (1 + pw[4]) * (c/2)
   
    for (i_duc=0; i_duc < 3; i_duc+=1)
        temp_duc_position_O[i_duc+1] = temp_duc_position_O[i_duc] + (1+pw[i_duc+4])*c
    endfor
   
        temp_duc_position_Y[1] = temp_duc_position_O[1] - (1 + pw[5]) * (c/2)
        temp_duc_position_Y[2] = temp_duc_position_O[2] - (1 + pw[6]) * (c/2)
        temp_duc_position_Y[3] = temp_duc_position_O[3] - (1 + pw[7]) * (c/2)

        temp_duc_position_Zr[1] = temp_duc_position_O[1] - (1 + pw[5]) * (c/2)
        temp_duc_position_Zr[2] = temp_duc_position_O[2] - (1 + pw[6]) * (c/2)
        temp_duc_position_Zr[3] = temp_duc_position_O[3] - (1 + pw[7]) * (c/2)
       
    Duplicate/O temp_duc_position_O temp_duc_position_O_reset, temp_duc_position_Y_reset, temp_duc_position_Zr_reset
    wavestats/Q duc_position_O
   
    temp_duc_position_O_reset = temp_duc_position_O - V_max
    temp_duc_position_Y_reset = temp_duc_position_Y - V_max
    temp_duc_position_Zr_reset = temp_duc_position_Zr - V_max

    variable i_F_duc
   
    wave FF_duc_Y, FF_duc_Zr, FF_duc_O, FF_duc 
    for (i_F_duc=0;i_F_duc<numpnts(temp_duc_position_O_reset);i_F_duc+=1)
        FF_duc_Y += atomic_frac_Y * ASF_Y* pw[i_F_duc] * exp(ci*xw*temp_duc_position_Y_reset[i_F_duc])*exp(-0.5*(xw*u_Y)^2)
        FF_duc_Zr += atomic_frac_Zr *ASF_Zr * pw[i_F_duc] * exp(ci*xw*temp_duc_position_Zr_reset[i_F_duc])*exp(-0.5*(xw*u_Zr)^2)
        FF_duc_O += 2*atomic_frac_O *ASF_O * 0.827 * exp(ci*xw*temp_duc_position_O_reset[i_F_duc])*exp(-0.5*(xw*u_O)^2)
    endfor
   
        FF_duc = FF_duc_Y + FF_duc_Zr + FF_duc_O


// ADSORBED ATOMS //
    wave FF_ads2, FF_ads
   
    FF_ads2 = ASF_Zr*pw[8]*exp(ci*xw*(pw[9]))*exp(-0.5*(xw*pw[10])^2)

        FF_ads = FF_ads2
       
// ADSORBED AND LAYERED WATER MOLECULES //
    wave FF_adw1, FF_adw2, FF_adw, FF_lw
    wave ASF_H2O

    FF_adw1 = ASF_H2O * pw[11] * exp(ci*xw*(pw[12]))*exp(-0.5*(xw*pw[13])^2)
    FF_adw2 = ASF_H2O * pw[14] * exp(ci*xw*(pw[15]))*exp(-0.5*(xw*pw[16])^2)
        FF_adw = FF_adw1 + FF_adw2
   
        FF_lw = ASF_H2O * exp(ci*xw*pw[17]) * exp(-0.5*(xw*pw[18])^2) / (1-(exp(-0.5*(xw*pw[19])^2)*exp(ci*xw*d_spacing_water)))

// MAGNIFICATION FACTORS //

    wave FFmag, FT_Cell, FB, FSpecular, Fsfact
    nvar cell_water_thickness, photon_wavelength, attenuation_length, r_e, suc_area
   
    FT_cell = exp((-8*pi*cell_water_thickness)/(xw*photon_wavelength*attenuation_length))
    FB = (1-pw[20])^2/(1+pw[20]^2-2*pw[20]*cos(xw*c))
    FSpecular = (4*pi*r_e/(xw*suc_area))^2

        FFmag = FT_Cell * FB * FSpecular * Fsfact


// FINALLY //
    wave FI_tot
   
        FI_tot = FFmag * magsqr(FF_uc*FF_ctr + FF_duc + FF_ads + FF_adw + FF_lw)
    return FI_tot
End
The problem with your fitting function is that you are trying to return a wave. Functions are expected to return one numeric value unless specified otherwise. If you need to, for some reason, specify the whole fit wave in one go, there is the option of an all-at-once fitting function. But also there you do not use return for the wave. You may want to take a look at the description of types and use of functions to clear these things up.
I recommend to at least go through the description of used-defined fitting functions by executing DisplayHelpTopic "Discussion of User-Defined Fit Function Formats"
Anyway, try to express your problem so that you get a defined output y for an input x, and return that. Or rewrite your function to an all-at-once type.