# How can I define the following function for arbitray values of the arguments?

Consider the following function

a[m_, l_] = 1/Log[m/l];


with this function I define

c[NN_, r_, l_] := cc /. Solve[NN == (1/a[1/r, l] - cc), cc][]


the things to keep in mind are these. The variable l is meant to be taken around 0.3 . NN is supposed to be an integer and for fixed l c will get smaller and smaller as NN-s increases. What I am interested at is to define a function that obtains for fixed l and r the value of NN that makes c smallest and positive. In order to do this I defined this

NNfun[r_, l_] := (NNN = 0; While[c[NNN, r, l] > 0, NNN++]; NNN - 1)


and it actually does the job well. Nonetheless I have a problem when I use this function to make a fit. I want to feed NNfun to another function in order to make a fit. The function I have in reality is more complicated than what I will show here, and the data I present here is not the data I have in but it is enough to make my point and to simplify the discussion. Thus, defining some random data points

data = Table[{RandomReal[], RandomReal[{-1, 0}]}, {i, 1, 30}];
dataError = Table[RandomReal[{0, 0.001}], {i, 1, 30}];


and defining the following function

fun[r_, l_] :=Sum[a[1/r, l]^n, {n, 1, NNfun[r,l]+1}]


I do a fit

NonlinearModelFit[data, fun[x, l] , {l}, x, Weights -> 1/dataError^2];
%["BestFitParameters"]
%%["EstimatedVariance"]


When I did this I obtained some numbers that made no sense to me. Eventually I think I discovered what is the issue behind. What I noticed is that

NNfun[r, l]


outputs -1 for an arbitrary value of r and l. The reason is that the While in the function does nothing for arbytray r and l so the output just takes the init value for NNN that I set which is 0 so the output is of course -1. as a consequence, when I evaluated

fun[r, l]


in the function I feed to NonlinearModelFit I actually am feeding 0, so that what I obtain is equivalent to doing

NonlinearModelFit[data, 0 , {l}, x, Weights -> 1/dataError^2];
%["BestFitParameters"]
%%["EstimatedVariance"]


So i finally arrive at my question, how can I define a function that just as NNfun above finds me the value of NN that makes c smallest and positive and that I can use in the fit function? that is that I can evaluate it for arbitrary r and l without obtaining nonsense. Or better, how can I modify my definitions so that I can use NNfun in the fit?

Here I wish to note that your data have non-positive values for the ordinate, while the function is positive in the given data range 0 < x < 1:
fun[r_Real, l_Real] := Sum[a[1/r, l]^n, {n, 1, NNfun[r, l] + 1}] 