module OW_mseq1
!		  Maximum length sequence assuming 2,3,5 distinct values
!
!       [ms]=MSEQ(baseVal,powerVal[,shift,whichSeq])
!
!       OUTPUT:
!       ms = generated maximum length sequence, of length basisVal^powerVal-1
!
!       INPUT:
!		  baseVal  -nuber of sequence levels (2,3, or 5 allowed)
!		  powerVal -power, so that sequence length is baseVal^powerVal-1
!		  shift    -cyclical shift of the sequence
!		  whichSeq -sequence istantiation to use 
!		  (numer of sequences varies with powreVal - see the code)

! (c) Giedrius T. Buracas, SNL-B, Salk Institute
! Register values are taken from: WDT Davies, System Identification
! for self-adaptive control. Wiley-Interscience, 1970
! When using mseq code for design of FMRI experiments, please, cite:
! G.T.Buracas & G.M.Boynton (2002) Efficient Design of Event-Related fMRI 
! Experiments Using M-sequences. NeuroImage, 16, 801-813.
 

subroutine owmseq1(baseVal,powerVal,shift,whichSeq,ms)

optional, integer				:: shift,whichSeq
integer 					:: baseVal,powerVal
integer					:: bitNum,i
integer, allocatable			:: register(:),tap(:),weights(:)
integer					:: ms(:)

if (not(present(whichSeq))) whichSeq=1
if (not(present(shift)))    shift=1

bitNum=baseVal^powerVal-1
allocate(register(powerVal)); register=1

if (baseVal==2) then
	select case (powerVal)
		case (2); tap(1).No=[1,2];
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/1,2/)
			end select	
		case (3)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/1,3/)
			case (2)
				allocate(tap(2))
				tap=(/2,3/)
			end select	
		case (4) 
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/1,4/)
			case (2)
				allocate(tap(2))
				tap=(/3,4/)
			end select	
		case (5)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/2,5/)
			case (2)
				allocate(tap(2))
				tap=(/3,5/)
			case (3)
				allocate(tap(4))
				tap=(/1,2,3,5/)
			case (4)
				allocate(tap(4))
				tap=(/2,3,4,5/)
			case (5)
				allocate(tap(4))
				tap=(/1,2,4,5/)
			case (6)
				allocate(tap(4))
				tap=(/1,3,4,5/)
			end select	
		case (6)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/1,6/)
			case (2)
				allocate(tap(2))
				tap=(/5,6/)
			case (3)
				allocate(tap(4))
				tap=(/1,2,5,6/)
			case (4)
				allocate(tap(4))
				tap=(/1,4,5,6/)
			case (5)
				allocate(tap(4))
				tap=(/1,3,4,6/)
			case (6)
				allocate(tap(4))
				tap=(/2,3,5,6/)
			end select
		case (7)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/1,7/)
			case (2)
				allocate(tap(2))
				tap=(/6,7/)
			case (3)
				allocate(tap(2))
				tap=(/3,7/)
			case (4)
				allocate(tap(2))
				tap=(/4,7/)
			case (5)
				allocate(tap(4))
				tap=(/1,2,3,7/)
			case (6)
				allocate(tap(4))
				tap=(/4,5,6,7/)
			case (7)
				allocate(tap(4))
				tap=(/1,2,5,7/)
			case (8)
				allocate(tap(4))
				tap=(/2,5,6,7/)
			case (9)
				allocate(tap(4))
				tap=(/2,3,4,7/)
			case (10)
				allocate(tap(4))
				tap=(/3,4,5,7/)
			case (11)
				allocate(tap(4))
				tap=(/1,3,5,7/)
			case (12)
				allocate(tap(4))
				tap=(/2,4,6,7/)
			case (13)
				allocate(tap(4))
				tap=(/1,3,6,7/)
			case (14)
				allocate(tap(4))
				tap=(/1,4,6,7/)
			case (15)
				allocate(tap(6))
				tap=(/2,3,4,5,6,7/)
			case (16)
				allocate(tap(6))
				tap=(/1,2,3,4,5,7/)
			case (17)
				allocate(tap(6))
				tap=(/1,2,4,5,6,7/)
			case (18)
				allocate(tap(6))
				tap=(/1,2,3,5,6,7/)
			end select
		case (8);
			select case (whichSeq)
			case (1)
				allocate(tap(4))
				tap=(/1,2,7,8/)
			case (2)
				allocate(tap(4))
				tap=(/1,6,7,8/)
			case (3)
				allocate(tap(4))
				tap=(/1,3,5,8/)
			case (4)
				allocate(tap(4))
				tap=(/3,5,7,8/)
			case (5)
				allocate(tap(4))
				tap=(/2,3,4,8/)
			case (6)
				allocate(tap(4))
				tap=(/4,5,6,8/)
			case (7)
				allocate(tap(4))
				tap=(/2,3,5,8/)
			case (8)
				allocate(tap(4))
				tap=(/3,5,6,8/)
			case (9)
				allocate(tap(4))
				tap=(/2,3,6,8/)
			case (10)
				allocate(tap(4))
				tap=(/2,5,6,8/)
			case (11)
				allocate(tap(4))
				tap=(/2,3,7,8/)
			case (12)
				allocate(tap(4))
				tap=(/1,5,6,8/)
			case (13)
				allocate(tap(6))
				tap=(/1,2,3,4,6,8/)
			case (14)
				allocate(tap(6))
				tap=(/2,4,5,6,7,8/)
			case (15)
				allocate(tap(6))
				tap=(/1,2,3,6,7,8/)
			case (16)
				allocate(tap(6))
				tap=(/1,2,5,6,7,8/)
			end select
		case (9)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/4,9/)
			case (2)
				allocate(tap(2))
				tap=(/5,9/)
			case (3)
				allocate(tap(4))
				tap=(/3,4,6,9/)
			case (4)
				allocate(tap(4))
				tap=(/3,5,6,9/)
			case (5)
				allocate(tap(4))
				tap=(/4,5,8,9/)
			case (6)
				allocate(tap(4))
				tap=(/1,4,5,9/)
			case (7)
				allocate(tap(4))
				tap=(/1,4,8,9/)
			case (8)
				allocate(tap(4))
				tap=(/1,5,8,9/)
			case (9)
				allocate(tap(4))
				tap=(/2,3,5,9/)
			case (10)
				allocate(tap(4))
				tap=(/4,6,7,9/)
			case (11)
				allocate(tap(4))
				tap=(/5,6,8,9/)
			case (12)
				allocate(tap(4))
				tap=(/1,3,4,9/)
			case (13)
				allocate(tap(4))
				tap=(/2,7,8,9/)
			case (14)
				allocate(tap(4))
				tap=(/1,2,7,9/)
			case (15)
				allocate(tap(4))
				tap=(/2,4,7,9/)
			case (16)
				allocate(tap(4))
				tap=(/2,5,7,9/)
			case (17)
				allocate(tap(4))
				tap=(/2,4,8,9/)
			case (18)
				allocate(tap(4))
				tap=(/1,5,7,9/)
			case (19)
				allocate(tap(6))
				tap=(/1,2,4,5,6,9/)
			case (20)
				allocate(tap(6))
				tap=(/3,4,5,7,8,9/)
			case (21)
				allocate(tap(6))
				tap=(/1,3,4,6,7,9/)
			case (22)
				allocate(tap(6))
				tap=(/2,3,5,6,8,9/)
			case (23)
				allocate(tap(6))
				tap=(/3,5,6,7,8,9/)
			case (24)
				allocate(tap(6))
				tap=(/1,2,3,4,6,9/)
			case (25)
				allocate(tap(6))
				tap=(/1,5,6,7,8,9/)
			case (26)
				allocate(tap(6))
				tap=(/1,2,3,4,8,9/)
			case (27)
				allocate(tap(6))
				tap=(/1,2,3,7,8,9/)
			case (28)
				allocate(tap(6))
				tap=(/1,2,6,7,8,9/)
			case (29)
				allocate(tap(6))
				tap=(/1,3,5,6,8,9/)
			case (30)
				allocate(tap(6))
				tap=(/1,3,4,6,8,9/)
			case (31)
				allocate(tap(6))
				tap=(/1,2,3,5,6,9/)
			case (32)
				allocate(tap(6))
				tap=(/3,4,6,7,8,9/)
			case (33)
				allocate(tap(6))
				tap=(/2,3,6,7,8,9/)
			case (34)
				allocate(tap(6))
				tap=(/1,2,3,6,7,9/)
			case (35)
				allocate(tap(6))
				tap=(/1,4,5,6,8,9/)
			case (36)
				allocate(tap(6))
				tap=(/1,3,4,5,8,9/)
			case (37)
				allocate(tap(6))
				tap=(/1,3,6,7,8,9/)
			case (38)
				allocate(tap(6))
				tap=(/1,2,3,6,8,9/)
			case (39)
				allocate(tap(6))
				tap=(/2,3,4,5,6,9/)
			case (40)
				allocate(tap(6))
				tap=(/3,4,5,6,7,9/)
			case (41)
				allocate(tap(6))
				tap=(/2,4,6,7,8,9/)
			case (42)
				allocate(tap(6))
				tap=(/1,2,3,5,7,9/)
			case (43)
				allocate(tap(6))
				tap=(/2,3,4,5,7,9/)
			case (44)
				allocate(tap(6))
				tap=(/2,4,5,6,7,9/)
			case (45)
				allocate(tap(6))
				tap=(/1,2,4,5,7,9/)
			case (46)
				allocate(tap(6))
				tap=(/2,4,5,6,7,9/)
			case (47)
				allocate(tap(8))
				tap=(/1,3,4,5,6,7,8,9/)
			case (48)
				allocate(tap(8))
				tap=(/1,2,3,4,5,6,8,9/)
			end select
		case (10)
			select case (whichSeq)
			case (1)
				allocate(tap(2))
				tap=(/3,10/)
			case (2)
				allocate(tap(2))
				tap=(/7,10/)
			case (3)
				allocate(tap(4))
				tap=(/2,3,8,10/)
			case (4)
				allocate(tap(4))
				tap=(/2,7,8,10/)
			case (5)
				allocate(tap(4))
				tap=(/1,3,4,10/)
			case (6)
				allocate(tap(4))
				tap=(/6,7,9,10/)
			case (7)
				allocate(tap(4))
				tap=(/1,5,8,10/)
			case (8)
				allocate(tap(4))
				tap=(/2,5,9,10/)
			case (9)
				allocate(tap(4))
				tap=(/4,5,8,10/)
			case (10)
				allocate(tap(4))
				tap=(/2,5,6,10/)
			case (11)
				allocate(tap(4))
				tap=(/1,4,9,10/)
			case (12)
				allocate(tap(4))
				tap=(/1,6,9,10/)
			case (13)
				allocate(tap(4))
				tap=(/3,4,8,10/)
			case (14)
				allocate(tap(4))
				tap=(/2,6,7,10)
			case (15)
				allocate(tap(4))
				tap=(/2,3,5,10/)
			case (16)
				allocate(tap(4))
				tap=(/5,7,8,10/)
			case (17)
				allocate(tap(4))
				tap=(/1,2,5,10/)
			case (18)
				allocate(tap(4))
				tap=(/5,8,9,10/)
			case (19)
				allocate(tap(4))
				tap=(/2,4,9,10/)
			case (20)
				allocate(tap(4))
				tap=(/1,6,8,10/)
			case (21)
				allocate(tap(4))
				tap=(/3,7,9,10/)
			case (22)
				allocate(tap(4))
				tap=(/1,3,7,10/)
			case (23)
				allocate(tap(6))
				tap=(/1,2,3,5,6,10/)
			case (24)
				allocate(tap(6))
				tap=(/4,5,7,8,9,10/)
			case (25)
				allocate(tap(6))
				tap=(/2,3,6,8,9,10/)
			case (26)
				allocate(tap(6))
				tap=(/1,2,4,7,8,10/)
			case (27)
				allocate(tap(6))
				tap=(/1,5,6,8,9,10/)
			case (28)
				allocate(tap(6))
				tap=(/1,2,4,5,9,10/)
			case (29)
				allocate(tap(6))
				tap=(/2,5,6,7,8,10/)
			case (30)
				allocate(tap(6))
				tap=(/2,3,4,5,8,10/)
			case (31)
				allocate(tap(6))
				tap=(/2,4,6,8,9,10/)
			case (32)
				allocate(tap(6))
				tap=(/1,2,4,6,8,10/)
			case (33)
				allocate(tap(6))
				tap=(/1,2,3,7,8,10/)
			case (34)
				allocate(tap(6))
				tap=(/2,3,7,8,9,10/)
			case (35)
				allocate(tap(6))
				tap=(/3,4,5,8,9,10/)
			case (36)
				allocate(tap(6))
				tap=(/1,2,5,6,7,10/)
			case (37)
				allocate(tap(6))
				tap=(/1,4,6,7,9,10/)
			case (38)
				allocate(tap(6))
				tap=(/1,3,4,6,9,10/)
			case (39)
				allocate(tap(6))
				tap=(/1,2,6,8,9,10/)
			case (40)
				allocate(tap(6))
				tap=(/1,2,4,8,9,10/)
			case (41)
				allocate(tap(6))
				tap=(/1,4,7,8,9,10/)
			case (42)
				allocate(tap(6))
				tap=(/1,2,3,6,9,10/)
			case (43)
				allocate(tap(6))
				tap=(/1,2,6,7,8,10/)
			case (44)
				allocate(tap(6))
				tap=(/2,3,4,8,9,10/)
			case (45)
				allocate(tap(6))
				tap=(/1,2,4,6,7,10/)
			case (46)
				allocate(tap(6))
				tap=(/3,4,6,8,9,10/)
			case (47)
				allocate(tap(6))
				tap=(/2,4,5,7,9,10/)
			case (48)
				allocate(tap(6))
				tap=(/1,3,5,6,8,10/)
			case (49)
				allocate(tap(6))
				tap=(/3,4,5,6,9,10/)
			case (50)
				allocate(tap(6))
				tap=(/1,4,5,6,7,10/)
			case (51)
				allocate(tap(8))
				tap=(/1,3,4,5,6,7,8,10/)
			case (52)
				allocate(tap(8))
				tap=(/2,3,4,5,6,7,9,10/)
			case (53)
				allocate(tap(8))
				tap=(/3,4,5,6,7,8,9,10/)
			case (54)
				allocate(tap(8))
				tap=(/1,2,3,4,5,6,7,10/)
			case (55)
				allocate(tap(8))
				tap=(/1,2,3,4,5,6,9,10/)
			case (56)
				allocate(tap(8))
				tap=(/1,4,5,6,7,8,9,10/)
			case (57)
				allocate(tap(8))
				tap=(/2,3,4,5,6,8,9,10/)
			case (58)
				allocate(tap(8))
				tap=(/1,2,4,5,6,7,8,10/)
			case (59)
				allocate(tap(8))
				tap=(/1,2,3,4,6,7,9,10/)
			case (60)
				allocate(tap(8))
				tap=(/1,3,4,6,7,8,9,10/)
			end select
		case (11)
			allocate(tap(2))
			tap=(/9,11/)
		case (12)
			allocate(tap(2))
			tap=(/6,8,11,12/)
		case (13)
			allocate(tap(2))
			tap=(/9,10,12,13/)
		case (14)
			allocate(tap(2))
			tap=(/4,8,13,14/)
		case (15)
			allocate(tap(2))
			tap=(/14,15/)
		case (16)
			allocate(tap(2))
			tap=(/4,13,15,16/)
		case (17)
			allocate(tap(2))
			tap=(/14,17/)
		case (18)
			allocate(tap(2))
			tap=(/11,18/)
		case (19)
			allocate(tap(2))
			tap=(/14,17,18,19/)
		case (20)
			allocate(tap(2))
			tap=(/17,20/)
		case (21)
			allocate(tap(2))
			tap=(/19,21/)
		case (22)
			allocate(tap(2))
			tap=(/21,22/)
		case (23)
			allocate(tap(2))
			tap=(/18,23/)
		case (24)
			allocate(tap(2))
			tap=(/17,22,23,24/)
		case (25)
			allocate(tap(2))
			tap=(/22,25/)
		case (26)
			allocate(tap(2))
			tap=(/20,24,25,26/)
		case (27)
			allocate(tap(2))
			tap=(/22,25,26,27/)
		case (28)
			allocate(tap(2))
			tap=(/25,28/)
		case (29)
			allocate(tap(2))
			tap=(/27,29/)
		case (30)
			allocate(tap(2))
			tap=(/7,28,29,30/)
		case default call seperr ("M-sequence is not defined for powerVal > 30"))
	end select
end if

if (not.present(whichSeq)) then
	call owexec_time (rnd)
	seed(1)=int(rnd/10.)
	call owexec_time (rnd)
	seed(2)=int(rnd+2000.)
	call random_seed(put=seed)
	call random_number(rnd)
	whichSeq=ceiling(rnd*size(tap))
else
	if (whichSeq > size(tap) .or. whichSeq<1) then
		write(0,*) ' WARNING: wrapping arround!' 
		whichSeq=mod(whichSeq,size(tap))+1;
   	end if;
end if;

allocate(weights(powerVal)); weights=0

do i=1,length(tap)
	weights(tap(i))=1
end do

do i=1,bitNum
   	! calculating next digit with modulo powerVal arithmetic
   	!   register, (tap(1).No)
	ms(i)=mod(sum(weights*register)+baseVal,baseVal);
  	! updating the register
  
  	register=cshift(register, -1)
  	register(1)=ms(i)
end do

!ms=ms(:);
if (not.present(shift)) then
   	shift=mod(shift, size(ms))
   	ms=cshift(ms,-shift)
end if

if baseVal==2,     
  	ms=ms*2-1;
else
  	call seperr('Error: wrong baseVal!');
end if

end subroutine

end module