Resulting polynomial surface, determined by linear least squares
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=R8), | intent(in), | dimension(1:long1, 1:long2) | :: | tab_in |
surface dont on a fait une approximation qui lui est soustraite |
|
| real(kind=R8), | intent(out), | dimension(1:long1, 1:long2) | :: | tab_out |
tableau résultant : surface |
|
| integer(kind=I4), | intent(in) | :: | long1 |
taille x |
||
| integer(kind=I4), | intent(in) | :: | long2 |
taille y |
||
| integer(kind=I4), | intent(in) | :: | nvarx |
degré du polynôme en x |
||
| integer(kind=I4), | intent(in) | :: | nvary |
degré du polynôme en y |
||
| integer(kind=I4), | intent(in), | optional, | dimension(1:long1, 1:long2) | :: | imask |
masque |
| logical(kind=I4), | intent(in), | optional | :: | verif |
dump |
|
| logical(kind=I4), | intent(in), | optional | :: | multi_thread |
use multithread? |
subroutine least_squares_tcheby(tab_in, tab_out, long1, long2, nvarx, nvary, imask, verif, multi_thread) !! Resulting polynomial surface, determined by linear least squares implicit none integer(kind=I4), intent(in ) :: long1 !! *taille x* integer(kind=I4), intent(in ) :: long2 !! *taille y* integer(kind=I4), intent(in ) :: nvarx !! *degré du polynôme en x* integer(kind=I4), intent(in ) :: nvary !! *degré du polynôme en y* real (kind=R8), intent(in ), dimension(1:long1, 1:long2) :: tab_in !! *surface dont on a fait une approximation qui lui est soustraite* real (kind=R8), intent(out), dimension(1:long1, 1:long2) :: tab_out !! *tableau résultant : surface* integer(kind=I4), intent(in ), dimension(1:long1, 1:long2), optional :: imask !! *masque* logical(kind=I4), intent(in ) , optional :: verif !! *dump* logical(kind=I4), intent(in ) , optional :: multi_thread !! *use multithread?* real(kind=R8), allocatable, dimension(:) :: var1, vec_x1, vec_x2, vec_hij real(kind=R8), allocatable, dimension(:,:) :: tab_tche1, tab_tche2, Jf real(kind=R8), dimension(1:(nvarx +1)*(nvary +1)) :: coeff_mxy integer (kind=I4) :: i, iu, nbvar, nbpts character(len=128) :: string if (nvarx==0 .and. nvary==0) then tab_out(1:long1, 1:long2) = 0._R8 return endif nbvar = (nvarx +1)*(nvary +1) if (present(imask)) then nbpts = sum(imask) else nbpts = long1 * long2 endif allocate( var1(1:nbvar) ) allocate( vec_x1(1:long1), vec_x2(1:long2) ) var1 = 0 do i = 1, long1 vec_x1(i) = -1. + (i - 1) * 2. / (long1 - 1)!real(i-long1/2, kind=R8)/(long1/2) enddo do i = 1, long2 vec_x2(i) = -1. + (i - 1) * 2. / (long2 - 1)!real(i-long2/2, kind=R8)/(long2/2) enddo call tab_tcheby(deg = nvarx, nx = long1, vec_x = vec_x1(1:long1), tab_tche = tab_tche1) call tab_tcheby(deg = nvary, nx = long2, vec_x = vec_x2(1:long2), tab_tche = tab_tche2) call tab_Jf_tcheby(nx1 = long1, & ! nx2 = long2, & ! nb_pts = nbpts, & ! nvarx = nvarx, & ! nvary = nvary, & ! nb_var = nbvar, & ! tab_tche1 = tab_tche1(1:long1, 1:nvarx+1), & ! tab_tche2 = tab_tche2(1:long2, 1:nvary+1), & ! tab_Jf = Jf, & ! imask = imask(1:long1, 1:long2), & ! multi_thread = multi_thread) ! allocate(vec_hij(1:nbpts)) if (present(imask)) then vec_hij(1:nbpts) = pack(tab_in(1:long1, 1:long2), mask = (imask(1:long1, 1:long2)==1)) else vec_hij(1:nbpts) = reshape(tab_in(1:long1, 1:long2), shape = [nbpts]) endif call moindres_carres_lineaire(nb_var = nbvar, & ! nb_pts = nbpts, & ! hij = vec_hij(1:nbpts), & ! beta = var1(1:nbvar), & ! Jf = Jf(1:nbpts, 1:nbvar)) ! deallocate(vec_hij) if (present(verif) .and. verif) then call coeff_poly_tcheby_xy_vers_poly_monome(var = var1(1:nbvar), & ! coeff_m = coeff_mxy(1:nbvar), & ! deg_x = nvarx, & ! deg_y = nvary) ! write(string, '(a,i3.3,a)') '(', nbvar,'E18.8,a)' call get_unit(iu) open(iu, file = 'verif_tcheby_vers_monome.txt') write(iu, trim(string)) (coeff_mxy(i), i=1, nbvar), " coeff tchebychev" ! à ce stade on a les coefficients var1(:) tels que : ! surface_approchante(x,y) = var1(1).t0(x)t0(y) +var1(2).t1(x)t0(y) +var1(3).t2(x)t0(y) +...+var1(nvarx+1 +1).t0(x)t1(y) +...+var1(nbvar).t_{nvarx+1}(x)t_{nvary+1}(y) write(iu,'(a)') 'pour vérifier avec gwyddion, il faut préalablement décocher "surface carrée", donner à l''image un facteur d''échelle x et y identiques, '// & ! 'et mettre un décalage pour que x soit compris entre -1 et 1 tout comme y. '// & ! 'coeff_mxy donne les coeffs de 1 x x^2 x^3 ... x^i y xy x^2.y ...x^i.y^j' close(iu) endif call tab_poly_tcheby(nx1 = long1, & ! nx2 = long2, & ! nvarx = nvarx, & ! nvary = nvary, & ! nb_var = nbvar, & ! tab_tche1 = tab_tche1(1:long1, 1:nvarx+1), & ! tab_tche2 = tab_tche2(1:long2, 1:nvary+1), & ! var = var1(1:nbvar), & ! tab_poly_tche = tab_out(1:long1, 1:long2), & ! multi_thread = multi_thread) ! deallocate( Jf ) deallocate( vec_x1, vec_x2 ) deallocate( var1 ) deallocate( tab_tche1, tab_tche2 ) return endsubroutine least_squares_tcheby