Yes, I had stripped out the Fortran from the geta subroutine downwards
because is wasn't very important in returning a test output matrix.
I think that I should make explicit what problem I believe I am seeing in
what I have coded, and what form of the output matrix I would prefer.
When I run the R code
> #Source positions in lens plane under annual parallax
> xypos_parallax <- function(year,ra,dec,ti,model_par) {
> if (!is.loaded('xypos_parallax_r')){
> dyn.load('parallax.so')}
> returned_data = .Fortran('xypos_parallax_r',
> as.integer(length(ti)),
> as.integer(year),
> as.double(ra),
> as.double(dec),
> as.double(ti),
> as.double(model_par$t0),
> as.double(model_par$tE),
> as.double(model_par$alpha),
> as.double(model_par$u0),
> as.double(model_par$piee),
> as.double(model_par$pien),
> matrix(0,nrow=length(ti),ncol=2) )[[12]]
> }
>
> model_par<-list(t0=5781.49344,
> tE=63.0814,
> alpha=1.2797,
> u0=0.0555,
> rho=0.0063,
> d=0.7119,
> q=0.0020496,
> piee=-0.900,
> pien=-0.800)
>
> ti<-seq(from=5770.0,to=5797.0,length=10)
> ra<-264.5590833
> dec<--27.13613889
> year<-2011
> y<-xypos_parallax(year,ra,dec,ti,model_par)
calling the Fortran subroutines;
> !**************************************************************************************************
> !XYPOS_PARALLAX_R: R wrapper subroutine
> !**************************************************************************************************
> subroutine xypos_parallax_r(k,year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,y)
> implicit none
> integer::year,i,k
> double precision,dimension(k)::ti
> double precision::ra,dec
> double precision::t0,tE,alpha,u0
> double precision::piee,pien
> double precision::y1,y2
> double precision,dimension(2,k)::y
> do i=1,k
> call xypos_parallax(year,ra,dec,ti(i),t0,tE,alpha,u0,piee,pien,y1,y2)
> y(1,i) = y1
> y(2,i) = y2
> end do
> return
> end
> !**************************************************************************************************
> !XYPOS_PARALLAX: Source positions under annual parallax
> !**************************************************************************************************
> subroutine xypos_parallax(year,ra,dec,ti,t0,tE,alpha,u0,piee,pien,xx,yy)
> implicit none
>
> double precision::ra,dec
> double precision::ti,t0,tE,alpha,u0
> double precision::piee,pien
> double precision::xx,yy
> integer::year
>
> double precision::qn, qe, thjd, t0dum
> double precision::dtau,du0
>
> thjd = ti
> t0dum = t0
>
> IF((piee .ne. 0) .or. (pien .ne. 0))then
> call geta(qn,qe,thjd,ra,dec,t0dum,year)
> !the following equations are still to be checked.
> dtau = piee*qn + pien*qe
> du0 = -piee*qe + pien*qn
> du0 = -du0
> ELSE
> dtau=0
> du0=0
> END IF
>
> xx = ((ti-t0)/te+dtau)*cos(alpha)-(u0+du0)*sin(alpha)
> yy = ((ti-t0)/te+dtau)*sin(alpha)+(u0+du0)*cos(alpha)
>
> return
> end
The matrix y which R returns is;
> -
*
> 0.09712576
*
> -0.036362910
*
> -0.17291343
*
> 0.067752761
> -0.08716623 -0.020590942
> -0.12098203 0.109764867
> -0.07617879 -0.003468398
> -0.07086033 0.149555601
> -0.06409356 0.015060440
> -0.02263996 0.187057573
> -0.05084316 0.035047997
> 0.02359192 0.222208831
But from Fortan, via
> do i=1,npoint
> write(*,*) y(1,i),y(2,i)
> end do
The following matrix is reported;
>
*
> -9.71257556964655877E-002 -0.17291342658178052
*
>
> -8.71662283852910891E-002 -0.12098202580263227
> -7.61787919231983052E-002 -7.08603339220746364E-002
> -6.40935606206668174E-002 -2.26399573881042940E-002
> -5.08431642024918945E-002 2.35919215412063593E-002
> -3.63629100550238102E-002 6.77527605328501203E-002
> -2.05909424429263285E-002 0.10976486692565257
> -3.46839842574317436E-003 0.14955560107302845
> 1.50604398227174099E-002 0.18705757282679777
> 3.50479966174077409E-002 0.22220883066716274
The second is the form of the matrix that should be stored in the returned
matrix y in the R function call, if possible.
Do you see how the elements are being stored in the wrong order in the
matrix that R is storing, versus the one that is read out in the Fortran? Is
it possible to get the form of the second matrix?
Thanks in advance.
--
View this message in context:
http://r.789695.n4.nabble.com/Trouble-returning-2D-array-into-R-from-Fortran-tp4646862p4647106.html
Sent from the R help mailing list archive at Nabble.com.
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
