function [free_physalia,mean_free,certainty_free_fract,encounter_risk]=PhysaliaSIM_mean_free_fun(minlon,maxlon,minlat,maxlat,data_longitude_free,data_latitude_free)
  
 
[xx,losdias,losanyos]=size(data_longitude_free);

for k=1:losanyos
for j=1:losdias
clear a, a=find(data_longitude_free(:,j,k)>=minlon&data_longitude_free(:,j,k)<=maxlon);
clear b, b=find(data_latitude_free(a,j,k)>=minlat&data_latitude_free(a,j,k)<=maxlat);
free_physalia(k,j)=length(b);


 % if length(b)>0,figure,medcoast,plot(data_longitude_free(a,j),data_latitude_free(a,j),'r.'),
 %   title(["day = ",int2str(j),", num colonies=",int2str(length(b))],'fontsize',20)
%"pause",pause,end
end %for j
 %"return",return

end %end k


%compute the mean numeric value (for all years within the database) of arrivals to our ROI
if losanyos>1
mean_free=nanmean(free_physalia); %nanmean lo calcula por columnas, no loop needed here
else
mean_free=free_physalia; %solo hay 1
end

%let's compute a certainty level [0 1] for the daily arrival to our ROI
% nanmean
 clear data_certainty
 data_certainty=free_physalia;
 data_certainty(data_certainty>0)=1;
 if losanyos>1
certainty_free_fract=nansum(data_certainty)/losanyos;
else %o sea, losanyos=1
certainty_free_fract=data_certainty;
end

%now the 'encounter risk' (RER)

plotlimit=max(mean_free)/100;
certainty_free_fract(mean_free<plotlimit)=0; %Para evitar que con num.colonias cero 
%en el plot, se vea que hay presence probability nocero.

rel_pres=(mean_free.*10)./(max(mean_free));
encounter_risk=rel_pres.*certainty_free_fract;



