[R] Trouble returning 2D array into R from Fortran

[paulfjbrowne]

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.