Below letter came from a Turkish professor via madsci. I called her and she says she is ready for help.
I want to be independent from not yet commercial softwares and i want to make the otf psf calculating software also.

It is a very good question. The photograph lenses are simulated by a lens
optical system model, which is based upon a generalized imaging system
model and coordinate system. The significant properties of the imaging
system can be formed from a real physical optical system. The object
plane and optical system are assumed to be plane and rotationally
symmetric respectively. The image formation process in an optical system
with incoherent illumination can be described by a convolution operation

Iimage(x,y,wavelength)=Point Spread function(x,y,wavelength){Convolution}
Iideal(x,y,wavelength)

where
Iideal(x,y,wavelength)={pi.T.Lplane(x/m,y/m,wavelength).R
(x,y,wavelength)}/{4.(f/Aeff)^2}

Iideal is the ideal, geometric image irradiance distribution, unaberrated
m is the magnification of the optical system
T is the transmittance of the optics
Lplane is the Lambertian object plane radiance distribution
R is the relative illumination factor
f/Aeff=effectivef-Number=0.5/Numerical Aperture

The Point spread function is the image of a point object formed by the
optical system. It contains all optical aberrations and diffraction in
the optical system. Iideal gives an ideal optical system without
aberrations or diffraction.

In same time, Image formation can be viewed as a linear filtering process.

FT{Iimage(x,y,wavelength)}=FT{Point Spread function(x,y,wavelength)}.FT
{Iideal(x,y,wavelength)}

FT{Iimage(x,y,wavelength)}=OTF(fx,fy,wavelength).F T{Iideal
(x,y,wavelength)}

Iimage(x,y,wavelength)= FT^-1{OTF.FT{Iideal(x,y,wavelength)}
where FT is the Fourier Transform Operator

You say "But If I am not wrong we have to change psf for each pixel color
and position."

You are right PSF(point spread function) depends on position and each
wavelength in coherent illumination

That is mathematics of the process

After that, by thinking the optical system (lenses), the image formation
can be achieved by convolving each section with its associated PSF
followed by summation of the individual results to obtain the final image
irradiance distribution. Alternatively, the final image airradiance
distribution can be obtained by filtering the spatial frequency
components of each section with its associated OTF, followed by inverse
Fourier transforming the filtered section transforms and summing the
results to obtain the final image irradiance distribution. The optical
system data are generated by lens design program for a given lens with
specific object plane distance. Lastly, the filtering is realized to form
the final image without distortion. The process is achieved step by step.

I hope, the answer can help you.
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
Best ,

Mustafa Umut Sarac